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Most of the time on this project was spent on the trajectory planning 
problem. As was guessed in the proposal the construction is equivalent to 
the classical spline construction in the case that the system matrix is nilpo- 
tent. If the dimension of the system is n then the spline of degree 2 n — 1 is 
constructed. This gives a new approach to the construction of splines that is 
more efficient than the usual construction and at the same time allows the 
construction of a much larger class of splines. All known classes of splines 
are reconstructed using the approach of linear control theory. Currently one 
paper is essentially finished and another will be finished before the first of 
the year. As a numerical analysis tool control theory gives a very good tool 
for constructing splines. However, for the purposes of trajectory planning it 
is quite another story. 

Consider the following simple situation x = u(t) and y = x(t ) and suppose 
that we want to track a signal that is of the form of a simple step function, 0 
to the left of zero and 1 to the right. Denote the function by s(t). Suppose we 
could track the signal exactly. Then y(t) = s(t ), y = s and finally u(t ) = s 
the derivative of a delta function. We have been able to prove using control 
theoretic techniques that the spline construction tracks not only the function 
but its first two derivatives. Thus the control will approximate the derivative 
of a delta function. Even if we disallow discontinuous trajectories the problem 
remains. If the function is piecewise analytic and continuous the control will 
still track a derivative of a delta function or the delta function itself. The 
problem is that the spline construction is too good. Some improvement can 
be achieved by matching the initial states and terminal states but there is 
still a problem. 

We’re in the process of developing schemes that will relax the constraints 
that require that the system pass exactly thorough a sequence of points, but 
only require that the system pass within a window centered at the point. I 
feel that this will give a satisfactory control. The relation of this scheme to 
asymptotic model following is interesting. In the classical asymptotic theory 
a signal is given and the control is constructed in order to bring the system 
and control together at infinity. Here we are asking that the signal and 
the control be brought together at a sequence of points and at fixed times. 
After some thought it is clear that the scheme is going to have problems. 
If we consider the rather mundane situation of trying to follow another car 
very closely we know that in order to maintain a very tight margin very 
high accelerations and decelerations are necessary. The exact same thing is 
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happening here. I along with the students will develop a theory to relax the 
constraints that impose the conditions of exact matching at specified times. 
Acceleration of thirty G’s are probably not suitable for a passenger aircraft. 

I have enclose four documents which contain reports of work done under 
this grant. The work is proceeding at a very good pace and I feel that we 
have made major accoomplishments during this year. 
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Abstract 

In this work, the relationship between splines and the control theory has been an- 
alyzed. We show that spline functions can be constructed naturaly from the control 
theory. By establishing a framework based on control theory, we provide a simple 
and systematic way to construct splines. We have constructed the traditional spline 
functions including the polynomial splines and the classical exponential spline. We 
have also discovered some new spline functions such as trigonometric splines and the 
combination of polynomial, exponential and trigonometric splines. The method pro- 
posed in this paper is easy to implement. Some numerical experiments are performed 
to investigate properties of different spline approximations. 


1. Introduction. 

Spline functions are well known and are widely used for practical approximation of func- 
tions or more commonly for fitting smooth curves through preassigned points. Spline tech- 
niques have the advantage over most approximation and interpolation techniques in that 
they are computational feasible. Most of the published spline algorithms are for polynomial 
splines and the vast preponderance are for cubic splines. There is a small but excellent lit- 
erature on the so called exponential splines and there is an even smaller literature on splines 
with more or less arbitrary nodal functions, [9, 3]. 

In this paper we will present a common frame work for splines that includes polynomial 
splines of all orders and generalized exponential splines of all orders. This common frame 
work is based on ideas from linear control theory. Let’s recall some basic ideas from control 
theory. A linear control system is a differential equation 

J t x(t) = Ax(t) + Bu(t) 


1 



where x 6 R n , u E R m and the matrices A and B are constant matrices of compatible 
dimension. The vector x is the state of the system and the vector u is the control. The idea 
is that we can use the control u to steer the state from point to point in the state space R n . 

We can think of the first component of x as representing the position of the system and for 
appropriate A the second coordinate is the velocity, the third acceleration, etc. A common 
situation, for example in air traffic control, is to specific the position that the system must 
be in at a sequence of times. So in fact what we have is a set of points through which 
the system must traverse at specified times. One could fit these points with a spline curve 
and then ask for the control that would move the system along that trajectory. In fact this 
can be done but we will show that the control law can be developed from natural control 
theoretic principles that will move the system through the points at the desired times and the 
resulting curve will be piecewise analytic and will have 2n — 1 continuous derivatives, i.e. a 
generalized spline. With this framework we can construct a wide variety of spline functions. 

If the matrix A is nilpotent then the resulting construction is just that for polynomial splines. 

If the matrix is 2 x 2 and one eigenvalue is zero and the other is a nonzero real number then 
the spline is the usual exponential spline. In general the nodal functions are the coordinate 
functions of the matrix function e At . 

In this paper we give a unified treatment of all of the common one dimensional spline 
functions using simple ideas from control theory. It is coming to be understood that there is 
a large overlap between linear control theory and elementary numerical analysis. Eigenvalue 
methods are know to be closely related to the theory of the matrix Riccati equation [2], there 
are close relations between observability and quadrature techniques [8], system identification 
and Prony’s method are very similar [I] and now we see that the spline constructions and 
basic linear controllability are manifestations of the same phenomena. 

In Section 2 we review basic material from the theory of linear control systems that is 
needed for the development and give a condition that characterizes the optimal control law 
that generates the spline functions. In Section 3 we give the details of the construction 
of spline functions using control theory and in Section 4 we classify the possible classes of 
spline functions that arise from the control theoretic construction. Irr Section 5 we examine * % 
in detail some of the particular classes from Section 4 and finally in Section 6 we present a m 
series of numerical examples comparing the various classes. 
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Let us divide [0, T] into n subintervals as 

0 — to <C fi <C ■ • • < tn — 1 ^ ~ T, 

and define h k = t k - tk-\, the length of the rth subinterval. Our goal is to find a control 
law u € C m-2 [ 0, T] that drives the system (2.1) from i(0) = x° to x(T) = x T such that the 
observation function y(t ) satisfies the interpolation conditions 

y(tk) = a k , k = 1, l‘ \ (2.4) 

» 

Furthermore, u(i) minimizes the functional 

[ T u{s) 2 ds. (2.5) 

Jo 

Such a control is called an optimal control. 

Definition. The system (2.1) is called controllable if for any x° , x 7 " , and r > 0, there is 
a u(f) such that, 

x 7 ='x(t) =5 e Ar S° +* r e A{r ~ s) bu(s)ds. 

Jo 

Theorem 2.1 : The system (2.1) is controllable if and only if 

rank (b, Ab,' • • , j4 m-1 6) = m. (2-6) 
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For the special matrix A as in (2.2), it is easy to verify that 




0 

0 

... 0 

1 \ 

0 

0 

... 1 
• • 

* 

0 

1 

• . 

• • • * 

* 

W 

* 

... * 

* ) 


(2.7) 


and hence the condition (2.6) is satisfied. From Theorem 2.1, the system (2.1) is controllable. 

i » ».< i T 

Theorem 2.2 : The system (2.1) is controllable, if and only if the matrix / e 3 bb e 9 as 

j o 


is invertible . 


For the special matrix in (2.2), we then define 

Mft) = {f e~ A *bb r e~ AT *ds)~ 1 . (2.8) 

Theorem 2.3 : When the system (2.1) is controllable, a control that moves the system from 
x(t) = p L to x(t) = pR given by 


uft) = Pe- AT \f e-^bPe-^'dsy'ie-^pR - e~ At -p L ), (2.9) 

minimizes the functional J{y ) = J i > 2 (s)<is among all controls that move the system from 
x(t) = p L to x(t) = p R . 

Theorem 2.4 : When the system (2.1) is controllable, a control u £ C m-2 [0,T] that moves 
the system from «(0) = x° , passing through <F xftk) = oik, to x(T) = x 1 is given by 

H — 1 m 

uft) = £ Pkfkft) + £ 7 .^( 0 . ( 2 ‘ 10 ) 

k = 1 1=1 


with 


f r t<t h 

t>t k , A: = 1, • • • ,n — 1, 

gi (t) = ^e A ^~% i = 1, • • • ,m, 


where e[ = (1, 0, • • • , 0), • • • , e£ = (0, • • • , 0, 1), and 0 k % 7 ; ’ are determined byn - 1 interpo- 
lation conditions (Fifth) = and m boundary conditions x(T) = iiF . Moreover , the control 
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(2.10) minimizes the functional J(v ) = / v 2 (s)ds among all functions v £ C fm “ 2 [0,T] that 
drives the system (2.1) from x(0) = x°, passing through irz(tk) = a to x(T) = x 1 . 

The construction of the optimal control is based on Hilbert space techniques and is based 
on writing the interior constraints in terms of a linear variety defined in terms of the functions 
fk(t) and the terminal constraints in terms of the functions Once the constraints are 

written in terms of linear varieties the form of the optimal control is clear based on the 
orthogonal complement of the intersection of the vaxieties. this is a standard technique and 
is found for example in [7] or [6]. The proof of Theorem 2.4 is based on two facts: (i) u(t) 
defined by (2.10) is m — 2 times continuously differentiable; (ii) /*’ s and g^s are n + m — 1 
linearly independent functions. In the following, we verify (i) and (ii). 

(i) From the construction (2.10), we only need to show that 

fk ) ( t k) = 0, r = 0, • • • ,m — 2, k = 1, •••,»» — 1. 

Indeed, for all k — 1, • • • , n — 1, r = 0, • • • , m — 2, 

lim fl r) (t)= lim eJ{-A) r e A ^- t) b=(-l) r e^A r b = 0, 

t-+t k - 0 K v 7 t— 7 7 

by virtue of (2.7). Hence fk(t) is m — 2 times continuously differentiable, and so is u(t). 

m 

(ii) Set = 0 for t 6 [0,T], i.e., 

1 = 1 

m 

F(t) = £ 7i ef e ^"-‘>6 = ^e A ^b = 0, 

1=1 

m 

with p = E Then we have F^(f) = 0 on [0, T], especially 

t=i 

F {r \t k ) = ?{-A) T b = 0, r = 0, • • • m — 1. 

Therefore, 

p{b,Ab,---,A m - l b) = 0 r . 

In light of (2.7), ~p is a zero vector ^ind consequently, 7,- = 0, t = 1, • • • ,m. So gfs are m 
linearly independent functions. 

Next we set 

/Wn-i(0 + Et Mt) = 0. [0,T]. (2.11) 

t'=l 
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If #,-1 / 0, then 


( 2 . 12 ) 


* I 

f n -l{t) = -~ £7 &(*)• 

Pn-1 i=1 

m t 

By the definition, f n -i{t) = 0 on So £7 .^<(0 = 0 for < e (i„-i,*n) which yields 

1=1 

7 - = o, i = 1, • • • ,m, and hence f n -x{t) = 0 by (2.12). This is a contradiction. Then (2.11) 
yields /3 n -i = 0, and consequently 7t - = 0, i = 1, • • * , m. So / n - 1 and <7t’s are m + 1 linearly 
independent functions. 

Continue the above procedure by adding f k s one by one, we are able to show that fk s 
and gf s are m + n — 1 linearly independent functions. 

3. Construction of splines by the control theory. 

Theorem 2.4 implies that an optimal control for the system (2.1) (with A given by (2.2)) 
is unique. But in general, /3t’s and 7 ,*’s in (2.10) are difficult to find, we then introduce 
a practical procedure to construct a control law that satisfies all the requirements. This 
control law actually leads us to a construction of spline functions. 

By the existence of a control law, there exists a set of points x 1 ,-*-,;? 1 1 with x x = 
a*, Jfc = l,---,n - 1 such that the solution of the system (2.1) satisfies x(t k ) = x k , k = 
0, 1, • * * , ti — 1, n. Bv virtue of Theorem 2.3, a control law that satisfies all the requirement 

can be defined piecewise as 

k = b • • • > ( 3 - 1 ) 

where u k (t) is given by (2.9) with t = t k -u t — t k) pL = PR = x * Then equations 

to find (n — l)(m — 1) unknowns in x 1 , ■ • * , x n ~ l (recall that x* = oc k k = l,***,n — 1, are 
known) come from ( n — l)(m — 1) continuity conditions on u(i), i.e., 

4 r) (^fc) = 4 r +i(^)> r = 0, • • • , m — 2, * = l,---,n-l. (3.2) 

From (2.9), 

u k (t k ) = (A T b) T e~ ATtk {f tk e~ Aa bb r e~ AT3 ds)~ 1 (e~ Atk x k — (3.3) 

= (A'b) T c- AT ‘'(f ,k ’\-^e- AT -ds)-\e- M ‘«?+' -e- A, ‘?y (3.4) 

• Next, we shall simplify (3.3) and (3.4). Toward this end, we introduce a change of variable 
s = tfc-i + s' into 

(jf* e ->u 6P" e~ AT, ds)~ l = (e'**- 1 e~ As ‘ We~ A J*‘ ’ ds t e~ ATtk ~ l )~ 1 
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k ~ x M(h k )e 


(3.5) 


= e 


A T t 


) 


where M(h k ) is defined by (2.8). Substituting (3.5) into (3.3), we have 

u[ r \t k ) = (A r b) T e~ AThk M(hk)(e~ Ahk x k - x fc ' 1 ). (3.6) 

Similarly, 

4+i(tJk) = ( A r b) T M(h k+1 )(e~ Ahk +'x k+1 - x k ). (3.7) 

Substituting (3.6) and (3.7) into (3.2) yields a linear system for (n - l)(m - 1) unknowns in 

-( A r b) T e~ AThk M(h k )x k ~ 1 + (A r b) T [e~ AThk M{h k )e~ Ahk + M(h k+1 )]x k 
-(A r b) T M(h k+1 )e- Ahk +'x k+l = 0, r = 0, ■ • • , m — 2, k = 1, (3.8) 


By virtue of the existence and uniqueness of the optimal control, the linear system (3.8) has 
a unique solution and hence its coefficient matrix is invertible. 

In order to solve (3.8), The following quantities needs to be calculated, A T , c~ Ah ( e~ A ht ), 
and M{h). Sometimes it is easier to use the Jordan matrix of A, denoted by A. There exists 
an invertible matrix Q such that A = and hence 

A r = QA T Q~\ e~ Ah = Qe~ Ah Q~\ e~ ATh = Q- T e~ ATh Q T . (3.9) 

M{h) = {£ Qe- A3 Q- i bb r Q- T e~ ATa Q T ds)- 1 = Q~ T M{h)Q~\ (3.10) 

where 

M(h ) = ( j\- Aa Q- l b{Q- l b) T e~ KTa ds ) _1 . (3.11) 

Substituting (3.9) and (3.10) into (3.8), we then have 

-(A r Q- x b) T e- KTkk M{h k )Q- l &- 1 + (A r Q- l b) T [e~ AThk M{h k )e- Ahk + M(h k+l )]Q~ 1 x k 
-{\ r Q- l b) T M{h k+l )e- Ahk +'Q~ 1 x k+l =0, r = 0, • ■ • , m — 2, k= 1, • • ■ , n - l.(3.12) 


Solving (3.8) or (3.12) for (n — l)(m — 1) unknowns in r 1 , •••,£" we then have the control 
u(<) defined piecewise by (3.1). The solution of the system (2.1) is thus given by 


f(0 = 


e At x° + 


f e A ^ 3 ^bu(s)ds 
Jo 


= Q(e At Q~ 1 x° + f e A ^ t ~^Q~ 1 bu(s)ds). 
Jo 


(3.13) 
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Note that x((t) = x, + i(t), i — l,***,m — 1. So the continuity of x((t) is continuity of 
x l+ i(£) for i < m. Further, continuity of Xi m+r ^(t) is continuity of r = 0, • • • ,m — 2. 

Therefore, the observation function y(t) = <?■ x{t) = xi(t) is a 2m — 2 times continuously 
differentiable function that satisfies the boundary conditions 

y (r) (0) =x° +1 , y^(T) = x? +1 , r = 0,---,m-l (3.14) 

and the interpolation conditions 

y{t k ) = x\, k = 1, • • • , n — 1. (3.15) 

Hence y(t) is a spline function. We see that from the control theory, we can derive quite 
general spline functions. Summing up, we have proved 

Theorem 3.1 : (1) There exists a unique function y(t) 6 C ,m “ 2 [0,T] that satisfies the 
boundary conditions (3.14) and the interpolation conditions (3.15); (2) y(t) is the first com- 
ponent of the vector function x(t) given by (3.13) in which u(s) is defined piecewise on each 
subinterval [<*- 1 , t *], k = 1, • • • , n, by 

u{ r) (t k ) = b r e- AT[t - t ^ ) M{h k ){e~ Ahk x k -x k ~ l ) 

= (Q~ 1 b) T e~ fLT ^ t ~ tk ~ l ^M(h k )(e~ Ahk Q~ l x k — <p -1 x fc-1 ), (3.16) 

where x k , k = 1, • • ■ , n — 1 are determined by solving the linear systems (3.8) or (3.12) (Note 
that x° , x" and x k , k = 1, ••*,«— 1 are given by the boundary conditions (3.14) an d the 
interpolation conditions (3.15)). 

In the next section, we will see that these splines can be piecewise polynomials, trigono- 
metric functions, exponentials or any combination. As special cases, we are able to recover 
classical polynomial splines (odd order) and exponential splines by properly selecting pa- 
rameters ai, • • • , a m in (2.2) for the matrix A. 

4. Classification of splines. 

The type of the splines is determined by its nodal shape functions. From the control 
theory, we are able to construct the nodal shape functions of splines. 

In order to see the kind of interpolation functions in x(t), we only need to consider one 
subinterval. Without loss of generality, we use the first interval (<o»^i) = (0, A) where the 
solution of the system (2.1) is given by 

x(t) = e M x° + f e A ^bu(s)ds. (4.1) 

JO 
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From Theorem 2.3, 


u(t) = b r e~ ATt { t e- i4 *S5 r e" jir- ds)" 1 (c" i ** f 1 - x°). (4.2) 

Jo 

Substituting (4.2) into (4.1), we have 

x (t) = e At x°+ f e A ^We- ATs dsM{h)(e- Ah x l - x°) 

Jo 

= g e A£ [Q- 1 f 0 + M(t)- 1 M(A)(e- A/ ‘g- 1 x 1 -C?- l f 0 )], (4.3) 

where M(h ) and M(/i) axe defined by (2.8) and (3.11), respectively. 

Theorem 4.1 : Let A be given by (2.2), let (pi{t), ■ ■ • ,p m {t)) be the first row of the matrix 
e At [I - M{t)~ l M(h)\ = Qe M [I - M (t)’ 1 M (h)]Q~\ (4.4) 

and let (qi(t), • • ■ , q m {t)) be the first row of the matrix 

e At M{t)- l M(h)e~ Ah = Qe At M{t)~ x M{h)e~ Ah Q-\ (4.5) 

Then for r = 0, • • • , m — 1, 

p! r) ( 0) = S i<r+1 , p < i r) (h) = 0, t = l,...,m, (4.6) 

<?j r) (0) = 0, qj r) (h) = Sj,T+i, j = (4.7) 


Proof : From Theorem 2.1 and (2.7), the system (2.1) is controllable. By virtue of 
Theorem 3.3, a control that moves x(f) from x(0) = x° to x(h) = x 1 is given by (2.9) (with 
t = 0, t = h, pi, = x°, pa = x 1 ), and consequently 

y( 0 = e[x{t) 

= ef (e At x° + J e A ^~ a ^bu(s)ds) 

= e{e At [I - M{t)~'M{h)}? + ef e M M{t)~ l M{h)e~ Ah x l 

= (Pl(0) • • • ^Pm(t))x° + (gi(t), • • • .^mW)® 1 

mm 

• '"•* = YLx°ipfit) + Y^x)qfit). . «"■ » * . (4.8) 

• »=i * • 

Choose x° = ei x 1 = 0 in (4:8), and we have 

Pi V \i) = y (r) (0> r = — 1. 
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Therefore for r = 0, • • ■ , m — 1, 

P,- r) (0) = p (r) (0) = *i r) (0) = Zr+l(0) = x° +1 = £,>+1, 

P; r) (h) = y (r) (h) = x[ r) (h) = x r+ i(h) = xj +1 = 0, i = 1» • * • *m. 

Then we have proved (4.6). The proof of (4.7) is similar. ■ 


We call pi,qi nodal shape functions by the characteristics (4.6) and (4.7). From (4.4) and 
(4.5), We see that the nodal shape functions are linear combinations of function entries of 
matrices e At and e A1 M(i) _1 . In order to see the type of functions in the spline, we only need 
to examine the entries of these two matrices. 

In the following, we classify the spline functions derived from control theory. This clas- 
sification is based on the spectrum of the coefficient matrix A of the system (2.1) under 
different circumstances. We shall concentrate on the case m = 2. The reasons are: (1) The 
general situation for large m is very complicated and is difficult to describe precisely. (2) 
The case m — 2 has almost all features for the general case. (3) From the practical point of 
view, the case m = 2 is the most useful and important case. Let 


A = 




/?,7 G R 1 . 


The eigenvalues of A are Ai = 7 + >/ 7 2 + $, A 2 = 7 — Vl 2 + ft- 


1. 7 2 -j- 0 > 0. There are two distinct real eigenvalues. Then the Jordan matrix A of A, 
the transformation matrix Q and its inverse are given by 


Then 


: Alt 0 \ 

0 e Ajt ) ‘ 


e A ‘ = 


From (3.11) (by changing h to t), we have 


a. )( _ 0 ( - i ’ i) ( 


® _Ai *. 0 




1.. 


( .(1 - e-* A “V 2A X (e - ( Al+Aj )‘ - 1)/(A X + A 2 ) 

(A 2 - Ax) 2 V \e-^ +x ^ - l)/(Ax + A 2 ) (1 - e- 2A -‘)/2A 2 


and hence 


(e A «* - e _Al ‘)/2Ai (e~ Aj< - e A “)/( Aj + A 2 ) 


e "W " (A 2 -A l} 2 { (e- A ^ - e A2t )/(Ax + A 2 ) ( e A >‘ - e - A ’‘)/2A 2 


(4.9) 

(4.10) 


ds 


(4.11) 


(4.12) 
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l.a. 7 ^ 0 , ^ 0 . In this case Ai, A 2 , — A 1? — A 2 are all distinct. We then have the 

exponential spline with basis functions given by linear combinations of e Alt , e -Alt , e Ajt , e -A2 ‘. 

l.b. 7 = 0, P > 0. In this case A x = — A 2 = the basis functions in l.a. degenerate. 
However, by applying the following limits 


e A ** — e Al * 


lim 

Aj Ai -f* A 2 

lim 


e A,«( e -(A 1 +A,)«_ 1) ^ _^ Ai< 
Ai+a 2 — *-o Ai 4" ^2 


= lim 


c -A!*_ c a 3 * 


, , = Urn 

A 2 — Ai Aj 4" ^2 Ai-fA ^— *0 A\ 4" A 2 


(4.12) becomes 


mm-' - 1 ( ( eA “ - e ' A ")/ 2A ' - (eA " 'l (4 13 ) 

M{ ) ”4A?( -te- x " e -*.‘)/2A, ) ’ ' J) 

Hence we have the exponential spline with basis functions given by linear combinations of 
e v^ e -v^ieVfr 

l.c. p = 0, 7 f 0. In this case, \ x = 0 (if 7 < 0), or A 2 = 0 (if 7 > 0). Again the basis 
functions in l.a. degenerate. Assume that Ai = 0, then A 2 = A = 27 . From the limits 

( gAi t — g— Ai tA 

lim Al— — } - = t, lim e A,t = 1 , 


A 1-.0 


2A X 


Ai — >0 


we have 


e At =( l i ) , M(tY 


0 e A 


1 _ 1 ^ e At 1 \ (4 

~ P \ e' xt - 1 (1 - e" 2A ‘)/2 ) ' 1 ’ 

(4.15) 


M e ~ U ~ l ) 

M( ) “ A 3 ^ i-e A ‘ (e A ‘-e - A ‘)/2 J - 

Therefore we end up with the exponential spline with basis functions given by linear com- 
binations of 1 , t , e 27t , e -27t . Later we shall further show that this is the classical exponential 
spline [9]. 

2 . 7 2 + f3 < 0. There are two complex eigenvalues: Aj = 7 + ioj, A 2 = 7 — iu>, where 

uj = a/-7 2 - $. 

2. a. 7 ± 0, P < 0. Evaluating (4.10) and (4.12), we have the exponential-trigonometric 
spline with basis functions given by linear combinations of e'^sinujt, e^cos ut, e-^sinwt, 
e~' r ‘ coswf...* * 

2 .b. 7 = 0, P < 0. Again this is a degenerated case where Ai = A 2 = iu = 
Therefore (4.10) is now 


e A< = 


cos uit + i sin ut 

0 


0 

cos u>t + 


i sin u>t ) ‘ 


(4.16) 
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Taking the limit 7 — ► 0 in (4.12), we then have 


At ^_, —1 ( sin utjuj — t(cosu>t + i sinwt) A 

6 4u i 2 \ — f(cosu;f — isinu>f) sin ut/u> < J ' 


(4.17) 


Hence, we have the polynomial- trigonometric spline with basis functions given by linear 
combinations of sin >/— 0t, cos y/— ~j3t, t sin t cos %/ — 

3. 7 2 + 0 = 0. In this case Ai = A 2 = 7. 

3.a. 7^0. We have non- degenerated Jordan form in this case, 

Therefore, 


(4.13) 


e A ‘ = 


e' rl te* 
0 




)-(:;)• 


(4.19) 


*M-* - T ) ( ‘( 7 ) (1/7 ’ 1) ( J )" 5 
- jC«- 1/7 r).*- 


(4.20) 


After some more detailed manipulation (see the next section) we can show that this is the 
exponential spline with basis functions given by linear combinations of e 7t , te 11 , e 7t , te 7t , 
similar to the case l.b. 

3.b. 7 = 0. In this case, A itself is a Jordan matrix. We compute directly the following 
quantities: 

A = (°, «*-(j ;), (4.21) 

' 3 /3 -t 2 / 2 


w-jff; r)CH(i 1 

(:S •',«). 


(4.22) 


(4.23) 


Then we have the polynomial spline with basis functions given by linear combinations of 
1 In the next section, we shall further show that this is the well-known cubic spline 

[4]. 

From the above discussion, we see that we may encounter all kinds of splines by varying 
parameters /? and 7. Two general cases axe l.a. and 2. a. where we have full sized exponential 
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or exponential-trigonometric splines. Degeneration occurs when zero or multiple eigenvalues 
appear. The extremal is the case 3.b. when both eigenvalues are zero. It is this extremal 
case that draws most of the attention. This is evidenced by extensive investigation regarding 
the cubic spline in the literature. Case l.c. also has been investigated from a different point 
of view. But we can hardly find any work regarding the other cases (except l.c. and 3.b.) 
listed above. 

The situation for m > 2 is similar. Let Ai, • • • , A m be eigenvalues of A. When Ai , • • • , A m ; 
— Aj, • • • , — A m are all distinct, we have the exponential spline with the basis functions given 
by linear combinations of 


e 


Ait 

* 


e 


-A,t 




See case l.a. When complex eigenvalues appear, we get the exponential-trigonometric splines 
with basis functions e Afc£ sin ujt, e~ Xkt cos ut (see case 2. a.). If we have multiple eigenvalues, 
the terms like 

te Xi , isincj£, ie“ A£ cosui, i 2 e A< , 

will appear in basis functions (see cases l.b., 2.b. and 3. a.). Finally, zero eigenvalues will 
introduce polynomials into basis functions (see case Lc.) and the extremal situation is that 
all eigenvalues are zero in which case we recover polynomial splines of order 2m — 1 (see case 
3.b.). 


5. Examples of splines. 

In this section, we shall work out in detail some classes of splines. We shall explicitly 
construct the nodal shape functions and the linear system needed to solve for the unknown 
parameters. 

L Our first example is the case 3.b. which turns out to be the classical cubic spline. We 
first construct the nodal shape functions. From Theorem 4.1 we need only to calculate the 
first row of matrix (4.4) and the first row of matrix (4.5). Let t — h in (4.22), and we have 


M{h) = 

Thus from (4.21) and (4.23), we have 


/t 3 /3 -h 2 / 2 \ 

-A 2 / 2 h ) 


-l 


_ 12 / h h 2 / 2 \ 

“ h* V ft2 /2 h 3 / 3 ) • 


n At 


- e*'M(t)-'M(h ) = ( J J ) - ( o 1 ) ( J 2 /2 ‘ 't ) T' ( 


h A 2 /2 \ 
h 2 / 2 A 3 / 3 ) 
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JftKSNftL M I » 


ranguftinr 





^ 1 t 'j _ ^ -2 t 3 /h 3 + 3 t 2 /h 2 t(-t 2 /h 2 + 2 t/h) ^ 

.# 

^ —2t 3 /h 3 + 3t 2 /h 2 t(-t 2 /h 2 + 2 t/h) ^ 1 -h ^ 


-2 t 3 /h 3 + 3 t 2 /h 2 t(t 2 /h 2 - t/h) 


Hence 


Pl (i) = l+2(i) 3 -3(|) 5 , M‘) = *[1 + ( {? - 2(f)]; (5-1) 

?1 (f) = -2(i) 3 + 3(i) 2 , fc (0 = t(j)[(j) - 1]. (5.2) 

They axe precisely the nodal shape functions for the Hermit interpolation. 

Next we set in (3.8), r = 0, x k — ( a k ,/3 k ) T , and h k = h k +\ = h (by this, we are using 
equally spaced intervals). The equation (3.8) is now 

- P e- ATh M{h)x k ~ l + P'[e~ ATh M{h)e~ Ah + M(h)]x k - Pm^)*-**#* 1 = 0, (5.3) 

for & = 1, •••,« — 1. Substituting 

„-ATK M(h ,( 1 0 \ 12 / h h 2 / 2 \ _ 12 / 1 h/2 \ 

e M{ti) ~y_ h i ) h 4 ^ h i/ 2 h 3 /3 ) ~ h 3 \ -h/2 -h 2 / 6 ) ’ 

M{h)e Ah = [e A ^ -J/q ) » 

1 h/2 \ ( l -h\ _12( 1 -V2^ 

h 3 { -h/2 h 2 / 3 


e _>l h M(h)e 
into (5.3) yieMs, 


12 / 1 h/2 ) f 1 -h 

h 3 l — /i/2 —h 2 /6 I l 0 1 


(A 1) ( ak 

[ h 2, h’[ h 




flk-\ + 4/?^ + /3jt+i — — (cKjt+i — Q;<:-x), k — 1, • • • ,n — 1. 


In the matrix form (5.4) is 


4 

1 

0 • 

■ 0 

0 

1 

4 

1 • 

■ 0 

0 

0 

1 

4 • 

■ 0 

0 

0 

0 

0 • 

• 4 

1 

0 

0 

0 • 

• 1 

4 



a 2 — <*o 
a 3 - ai 
a 4 - a 2 

«n-l - On-3 

a„ - a n _2 



(5.5) 
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This is precisely the same linear system as we construct the cubic spline. After solving /3k' s 
s from (5.5), the desired cubic spline can be expressed piecewisely on the subinterval 
k = 1,* • • ,n as 

y(t) = — tfc-x) + /3k-iP2(t — tk-i ) + cikqi(t — 


where Pi,P 2 > 4 i >92 are given by (5.1) and (5.2). 

2. The second example is the case l.c. which is the classical exponential spline [9]. The 
nodal shape functions are again derived from the first rows of matrices (4.4) and (4.5). We 
use the Jordan form, let Aj = 0, and A 2 = A (4.9), and we have 


«-(;;)■ e-'-ic -/)• 


Using (4.14) we can calculate (by setting t = A), 


M(h) = A 3 
where 


Xt 


A t 


- 1 


-1 


r M -l (1 — e~ 2Xt )/2 


( 


(1 + e _A ")/2 1 

Xh(l - e 


A/i(l + e~ Xh ) — 2(1 — e~ xh ) 
Recall (4.14) and (4.15), and we have 


- C(X) 1 2 \un _ 


2A 3 


Qt Ki [I - M{h)\Q 


-1 


1 0 

0 e A ‘ 


[/- 


gw 

A 3 


Xt 

1 — e At (e 


= (:!)( 

( 

- (; 

( (l + e - AA )/2 1 WA -1 \ 

^ 1 Xh(l -e~ Xh )~ l ) V 0 1 / 

_ ^Pi(A;t) P2(A;0j 


e _A< — 1 \ . 

A ‘-e- A ‘)/2 ) ’ 


(l+e- AA )/2 1 U/A -l\ 

1 Xh(l — e _AA ) _1 y A \ 0 1 ) 

(e At - 1)/A \ C(A) f 1 1 W Xt e- At - 1 

0 e At ) A 4 V 0 At ) \ 1 -e A ‘ (e Al - e" A *)/2 

(1 + e- AA )/2 


where 


Pi(A; t) = 1 - 


Ai(l + e -AA ) - (1 - e" A *) - e A (‘" A ) + e 
Ah(l + e _AA ) — 2(1 — e~ xh ) 


-Xt 


(5.6) 
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Pj(A;0 = 


.At 


1 — e~ Xh t 1 — e xt 1 ~ (1-e-AT.j e At — 2 + e 


1 + e 


— Xh 




.-At 


A/» 


AA 


'(t 4 


AA 


-) 


AA 


AA 


(5.7) 


Qe M M{t)~ l M(h)e~ tih Q~ x 

C{ A) / 1 1 W 1 0 W A* e~ At - 1 \ 

~ A 4 V 0 A ) ^ 0 e Xt ) \ e~ Xt - 1 (1 - e~ 2A ‘)/2 ) ' 

( (l+e- AA )/2 1 Wl 0 W A -1 \ 

\ 1 AA(1 — e -AA ) -1 / \ 0 e~ xh j\0 1 ) 

C{ A) / 1 1 W At e _At — 1 \ ( (l+e- AA )/2 e" AA 

“ A 4 (o A j ( l-e A ‘ {e xt - e~ Xt )/2 ) { 1 AA(e A/l - l)" 1 



^ ?i(A;<) ?2(A;t) \ 


where 


?2(A;0 = 


?i(A;f) = 1 -pi(A;t), 


r l — e~ xh .t 1 — e'* e AAjj 

— ( l + — J -- 


a Xt 


Xh 


1 + e“ A/l - 2 


A h 


0 X t 


Xh 


— 2 -f* c 
~Ui 


-A t 


(5.8) 

-]• (5-9) 


(5.6) - (5.9) are nodal shape functions for the classical exponential spline. 

Next we set in (3.12), r = 0, x k = (a*,/?*) 7 *, and = Ajt+i = h . The equation (3.12) is 


now 


-{Q- x b) T e- KTh M{h)Q- l x k ~ x + {Q- x b) T [e- tJh M(h)e- Kh + M{h))Q~ 1 x k 
-{Q- x b) T M(h)e~ Kh Q- l x k + l = 0, k= 1,- (5.10) 


* Substituting 

(Q-'b) T = ±(-1,1), 


e -A/ *M(A) 


M(h)e~ xh 


= C( A) 

- • 

= C(A) 



(1 + e -AA )/2 
~-Xh 


= [e- M M(A)] T = C(A) ( 


(1 + e~ Xh )/2 1 

1 AA(1 - e -A/l ) -1 


Xh(e xh - I) -1 J ’ 

(1 + e~ Xh )/2 c~ Xh \ 

1 AA(e AA - I)" 1 ) ’ 


16 



M(h)e-^ = C(A) ( J e _° sl ) ( (!+'"“ 


)/2 


\h(e 


e~ Xh \ 
A^_i)-l J 


(1 + e~ Xh )/2 


= C-(A) 

into (5.10), canceling C r (A)/A 2 , we have 

+(-i,i) 

-(-i,i) 


(1 + e~ Xh )/2 


,-xh 


e -XK 

A ke~ Xh {e Xh 


- i)-' ) 1 


h, 




t. Xh(e xh 

l + e~ Xh l + e~ xh ) f A — 1 A / 

-AA _l „-2A/>n _ <,-A/i)-l I l o 1 11 


1 4- e A/ * Xh(l + e 2A/ *(1 — e~ 


<*k 

A 


(1 + e~ xh )/2 


,—Xh 


1 


A/i(e A " - 1) 


■')(-■)( 


a * +1 |=0 
ft+i j U ’ 


or 


1 _ e -™ _ 2Xhe- x \ AA(1 + e _2A/l ) _ Xh 

_ e _AA\ (^-! + &+W + l j _ e — A A ( J + e )]& 


2(1 


1 - e- AA 

~ A 7 (tAfc+l Orjt_l )i fc=l,*’",U 1. 


(5.11) 


This is the linear system for the exponential spline, it can be written in the matrix form as 


( a(\;h) 

b{\;h) 

0 

0 

0 ) 


0i \ 

b(\-,h) 

a(A; h) 

b(\;h) 

0 

0 


02 

0 

b( A; h) 

a(A; h) 

0 

0 


03 

0 

0 

0 

a(A; h) 

b( A; h) 


0n-2 

0 

0 

0 

... 6(A; A) 

a(X-h)) 


{ 0n- 1 / 


3 

h 


( 0.2 OtO \ 

(*3 — Ct i 

a 4 - a 2 

a„_i — ar n _ 3 

V <*n-2 / 


l b( A; /i)A> \ 
0 
0 


0 

V b(X;h)0 n J 


where 


jl\ fi AA( 1 + e~ 2A/l ) - (1 - e- 2A *) 1 - e- 2AA - 2AAe- A * 

a(A;A) = 6 — , b(X, h) = 3- 


(5.12) 


Xh( l-e~ Xh ) 2 ’ ~ AA(1 — e~ Xh ) 2 ’ 

Again, after finding 0^s, the spline can be expressed piecewise by the nodal shape functions. 
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It is interesting to examine the limiting case for the exponential spline obtained above. 
Applying the L’Hospital’s rule three times, we are able to verify that 


lim a(A; h) = 4, lim 6(A; A) = 1. 


(5.13) 


We then recover (5.5), the tridiagonal systems for the cubic spline. Further we can verify 
that (by successively using the L’Hospital’s rule) 

E 3*(* " f * - - 1 - - » 3 - < S > 3 = - 2 '|) 3 + 3 < S >’ = »<*>• 


qi(t) is one of the nodal shape functions for the cubic spline given by (5.2). The other three 
nodal shape functions can be verified similarly. So the cubic spline is the limiting case for 
the exponential spline 4.I.C. when A — ♦ 0. This is not a surprise from the following limit 
regarding matrix A of the system (2.1), 


lim 

A— »0 




(5.14) 


where the left hand side is the matrix A for the case l.c., and the right hand side is the 
matrix A for the case 3.b. We can also examine the limit when A — ► oo where 


lim a(A; h) = 6, lim 6(A; h) = 0, (5.15) 

A — *oo A— ►oo 

lim pi(A;t) = 1 - lim ? X (A; t) = -, (5.16) 

A— *>oo n A— *-00 (I 

lim p 2 (\',t) = 0 = lim q 2 (X;t). (5.17) 

A— ►oo A— ►co 

Substituting (5.15) into (5.12), we then have /?* = (a* + i —atk-i)/2h , k = 1, • • • , n — 1, which 
is the central difference scheme. We see that (5.16) gives the nodal shape functions for the 
linear interpolation. In order to verify (5.17), it is convenient to rewrite 


P2(A;<) = 


92(^5 0 «- 


,-A/i 


(Xh + 2)e 


-Xh 


(?- 


H(l-e" AA ) + 


,-Xh 


2A 3 


-2 


.A (t-h) 


o-At 


Xh-2 1 — e~ Xh (Xh — 2)[A/i(l + e _AA ) — 2(1 — e _M )] v A 1-e 
1 C(A) 

A 

C( A) 

2A 3 


ITT*) 


rwr(e- A ‘-2)], 


1-e 


H(1 - e~ Xh ) + 


1 _ g— At g-A/i _ gA(t— fc) /j( gA(t-A) _ 2e -AA + e -A (t+M) 


+ 


+ 


0 -Xh 


So the linear spline (the piecewise linear interpolation) is the limiting case for the exponential 

spline 4.1. c. when A — ► 00 . 
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3. The third example is the case 2.b. when 0 = — u> 2 , 7 = 0 and 



Denote 

Ci = u 2 h 2 — sin 2 cj/i, C2 = ^ sin 2uh + uih, C3 = sin 2 u>h, C4 = - sin 2uh — ujh. 

Following the same procedure as the first example, we have the nodal shape functions as 
following 


Cipi(t) = Ci cosut + Casino;* — utcosut) — Czutsinut, 

CiP2{t) = (jjh 2 sin u/t — C^t cosut + C^sinwt, 
qi{t) = 1 -Pi(t), 

C iq2 (t) = h sin uh(sinwt — u>t cosut) — t sinu>t(sinuh — uh cos u>h). 

When evenly spaced intervals are used, the tridiagonal system for unknowns 0 k is given by 


b(u] h)0k-i + a(w; h)0k + b(u\ h)0 k +\ = c(u;; h)(a k +i - ajfc-i), k = 1, • • • , n, (5.18) 


where 


a(u; h ) = 2 u>h — sin 2uh, 


b(ui ; h ) = sinu >h — uhsinuh, c(u>\h) = u 2 h sin ujh. 


It is easy to verify that, 


lim ¥r^ 

w-»o (u;/i ) 3 


= 4, 


36(cj; h) 

lim 

w-*0 (u/l) 3 


= 1 , 




3 

h' 


which are the coefficients for the tridiagonal system of the cubic spline. We can also verify 
that at the limit u> — ► 0, the nodal shape functions p,-(u;;t) and i = 1,2, have the 

relative nodal shape functions of the cubic spline as their limit when u — * 0. Indeed, we 
recover the cubic spline from this trigonometric spline when w-»0. 


4. The fourth example is the case 3. a. where 0 = — 7 2 and 



Denote 

C = 1 - 2e~ 2lh + e" 4 ^ - 47 2 h 2 e~ 2yh . 
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Computing the first rows of matrices e At — e At M(t ) 1 M(h) and e Ai M(t ) 1 Af(h)e Af> , we 
have the nodal shape functions 

Cpi(t) = e-^ 2 /l+t) ( 27 h - 2y 2 h t - 7 * + 2 ^ht - 1) + e~*(l + 7 i) 

+c -r(a *-*)( T t + 2t 2 ht - 2'yh + 2 7 2 h 2 - 1) + e -M^-0(i _ 7< ) > 

C7pj(<) = t(e- 7( ‘ ,A - f) + e" 7 ‘) + e~' l(2h ~ t \2 1 ht - 2 ~fh 2 - t) + e~^ 2h+t) {2 7 h 2 - < - 2 7 ft<), 

?i(0 = i-pi(<)> 

Cq 2 (t) = e-^-^ih - t - 2 7 ht) + (e~^ 3h+t) + e -^-0)( f _ h) + e" 7 ^ (2 7 hf + h - t). 


Evaluating (5.3) in the current case, we get the tridiagonal system (5.18) with 

a(u;; h) = a( 7 ; h) = e 7/l (l — e -47/ * — 4 7 he 27/l ), 
b(uj; h) = b( T h) = e- 27A (l+ 7 h)-(l - 7 h), 

c(w; /t) = c( 7 ; h) = 7 2 A( 1 - e -2Vl ). 


Again we can verify that the cubic spline is the limiting case for this exponential spline when 
7 — ► 0 . 

5. Our last example is a case for m = 3 when 


0 1 0 
A = | 0 0 1 
0 0 0 


The system (2.1) with A given by (5.19) produces the quintic spline. In this case, 


e At = 


1 t t 2 [2 \ 
0 1 t 
0 0 1 / 


/t 5 /( 20) —t 4 /8 t 3 /6 

MU)- 1 = I e-^We-^’ds = -t 4 / 8 t 3 / 3 -t 2 /2 

Jo V t 3 / 6 -t 2 / 2 t 


M{h) 


[M(h)- 1 ]- 1 = 3 


/6 -t 2 S 2 

240 /h 5 120 jh 4 20/A 3 \ 
120 /h 4 64/A 3 12/A 2 . 

20/ A 3 12/A 2 3/A / 


(5.19) 


(5.20) 

(5.21) 

(5.22) 


Using (5.20) - (5.22) to compute the first rows of matrices e At - e At M{t) l M(h ) and 
e At M(t)- 1 M(h)e~ Ah , we then have the nodal shape functions for the quintic interpolation: 

p,(() = ^/!V + 3*.f + 6( J ], „(() = ^[A J -3A( + 6( J ], 


h 5 ' 
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Pl(t) = 
ps (0 = 


( ^T~l ht + 3t>], 9 ! (i) = ^(Mi-*)-3( ! ], 


(h - t) 3 t 2 

2h 3 


$$(*) = 


t 3 {t - hf 

2h 3 


In order to compute the parameters for the optimal control, we set in (3.8) r — 0, 1 x k — 
(ajfc,/?fc, 7 jt) r , and = hk + i = h. Then except (5.3), we also have 


- ( Abfe- ATh M(h)x k - 1 + (A5) r [e-^'‘M(/ l )e-^ + M(^)]x fc - (Abf M(h)e~ Ah x k+1 = 0, 

(5.23) 

for k = 1 , •••,»» - 1. Substituting ( 5 . 20 ) - (5.22) into (5.3) and (5.23), after some tedious 
symbolic manipulation, we have the following linear system, 

20 

8(/3fc+i - 0k- 1 ) + h{~ 7*-i + 6 7; - 7fc+i) = y(A-i - 2/* 4- A+i); (5.24) 


7/3*:_i + 16/3fe + 7/3jt + i + h(~fk-\ - 7fc+i) = ~^~(A+i - A-i)» (5.25) 

for Jb = 1 , • • • ,n — 1. If we arrange the unknowns as (/?i, /» 7 i, * * • , A»-i> ^ 7 n-i)> we w dl get 
a linear system with a banded 6 -diagonal coefficient matrix; if we arrange the unknowns 
as (/?!,- ••,/?n-i,/i7i,---,/i7n-i) = (/ 3 T , /r? T ), we will have a linear system with a block 
tridiagonal coefficient matrix. 


where 


' s 

E T ' 

r 


r i 

/ 

. 8E 

H 

h~f 


. 9 . 


' 16 

7 

0 

... o 

0 ‘ 


* 6 

-1 

0 

0 

1 

o 

7 

16 

7 

... 0 

0 


-1 

6 

-1 ••• 

0 

0 

0 

7 

16 

... 0 

0 

, H = 

0 

-1 

6 ••• 

0 

0 

0 

0 

0 

••• 16 

. 7 


0 

0 

o ... 

6 

-1 

1 

o 

0 

0 

... 7 

16 _ 


0 

0 

0 ••• 

-1 

6 _ 


(5.26) 


' 0 

l 

0 •• 

0 

o 

-1 

0 

1 •• 

0 

0 

0 

-1 

0 •• 

0 

0 

0 

0 

0 •• 

0 

1 

t 

O 

0 

0 •• 

-1 

o 
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’ 15 (/ 2 - fo)/h - 10o - h~f 0 


20 (/o — 2/i + + 8/?o + ^7o 

15 (/ 3 -/i)//> 


20(/i — 2/2 + f 3 )/h 

15(/4-/ 2 )/A 

, 9 = 

20(/2 — 2/3 + fi)jh 

15(/„_i - f n . 3 )/h 


20(/„_3 — 2/ n _2 + fn-l)/h 

. 15 (/„ - fn-2)/h ~ 10 n + h^n _ 


. 20(/ n _2 — 2/„_x + fn)/h — 80n + h~f n _ 


In general, if A is a m x m nilpotent matrix with l’s on the super diagonal and 0’s 
elsewhere, we shall recover all odd degree polynomial splines (with degree 2m — 1). 

6. Numerical experiments. 

In this section, we test the behaviors of different splines numerically. Equally spaced 
intervals axe used for all computations. 

Example 1. Comparison of the cubic spline with the quintic spline. 

Test function 1. 

( 0 -1 < t < 0 

m =1/2 t = o 

[ 1 0 < t < 1 

For the cubic spline, we pose, in (5.5), the boundary conditions: 

Po = o = p n ; 

and for the quintic spline, we pose, in (5.25), the boundary condition: 

0o — 0 = 0m 70 = 0 = 7n. 

Recall that /?, and 7 j are the coefficients for the first and the second derivatives, respectively. 
The spline functions are then constructed for h = .2, h = .1, h — .05, and h = .025. Graphs 
are plotted in Figure 1(a), 1(b). We see that the qualitative behavior of the two splines are 
• ' almost same, but the quintic spline has a little better accuracy. 

One interesting phenomenon is that the mesh refinement does not effect the maximum 
overshoot of the spline approximation. Since this is very similar to the Gibbs phenomenon 
for the Fourier series, we term it as “Gibbs phenomenon” of splines. In fact, all spline 
functions have this property. 

Test function 2. 

g(t) = e~ wt \ -l<t< 0. 
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For the cubic spline, we pose the boundary conditions: 

/?0 = — 30e 10 , /?„ = 0; 

and for the quintic spline, we pose the boundary conditions: 

A> = — 30e 10 , 0 n = 0, 70 = 960e 10 , 7n = 0. 

We use the mesh size h = .2, and plot the graphs in Figure 1(c), 1(d). We see that the 
quintic spline gives much better approximation in the neighborhood of x = — 1 since it has 
the correct concavity information at x = — 1 which the cubic spline does not have. 

Example 2. Properties of the classical exponential spline case l.c. 

The test function is the same f(t) as in Example 1. We have observed that for small 
parameter A, the behavior is much like the cubic spline. This is not surprise from (5.14). The 
interesting fact is: For the moderate A, the graph is very much the same as the cubic spline 
(Figure 2(a)). If we fix the parameter A and refine the mesh, we observe Gibbs phenomenon 
as in the cubic and the quintic splines. But if we fix the mesh (here we choose h = .1) and 
increase the parameter A, we see that the approximation converges to the piecewise linear 
function (Figure 2(b), 2(c), 2(d)). This confirms our theoretical analysis made in Section 5. 

Example 3. Properties of the exponential spline case 3. a. 

We use the same test function f(t) as in the examples 1, 2. 

For small parameter 7, the approximating feature of this spline is also like the cubic spline 
including the Gibbs phenomenon. But when we fix the mesh (here h = .1) and increase the 
parameter 7, an unexpected wiggling appears at t = 0 (Figure 3(a)-3(d)). 

Example 4. Properties of the exponential spline case l.a. 

Here we choose 



and we have Aj = —1, A 2 = —2. Again, the testing function is f{t) as in the previous 
examples. 

We plot the approximation for h = .1, h — .05, h = .025 in Figure 4(a), 4(b), 4(c), 
respectively, again, we observe the similar behavior as that of the cubic spline. 

Conclusions 

1. Gibbs phenomenon exists for all splines. 

2. The quintic spline is recommended if the concavity is important. 
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3. From the approximation point of view, the classical exponential spline is preferred 
when the function has points of discontinuity. 

A final remark. 

For the discussion purpose, we constructed spline approximation in this paper by intro- 
ducing the nodal shape functions which is not necessary in practical computation. From the 
framework we have established based on the control theory in Section 3, all we need to do 
is: providing the matrix A, the vector b to the linear system (3.8), solving (3.8) numerically 
to obtain x k 's, and hence the control law u(i) (see (2.9)). After we have the control u(f), the 
expected spline function is given by the first component of i(t) defined by (3.13). Based on 
our analysis, we are able to choose different splines by simply selecting entries of the matrix 
A. 

The significance of this investigation is two fold: first, it exposes the relationship between 
two important fields - control theory and spline approximations. This enables us to discover 
new spline functions and to investigate, systematically, the properties of the spline approx- 
imations. Secondly, it provides a practical way to construct different splines from a same 
simple framework. From our experience, we feel that this construction is more natural and 
easier than the traditional approach. 
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A BRIEF SURVEY OF CONSTRAINED MECHANICS AND 0 , zf 

VARIATIONAL PROBLEMS IN TERMS OF DIFFERENTIAL FORMS / " ' 


Robert Hermann 

Ever since my graduate student days, I have been impressed an d 
influenced by the elegance and systematization of Mechanics and 
Variational Calculus contained in Elie Cartan's book "Lecons surles 
Invariants Integraux". In the period 1959-69, 1 expended considerable 
effort in the development of Cartan 's point of view in many books an d 
articles. In this paper (which will appear as a Chapter in "Interdisciphnay 
Mathematics ", v. 30), I will give a quick development of some of the 
material in my books " Differential Geometry and the Calculus of 
Variations" and "Geometry, Physics and Systems". 


Another purpose in developing this geometric form of the Equations 
of) Mechanics in this Volume is that it fits in with my strategy of 
investigating mechancis with 'singular' features, such as Delta Functions, 
Discontinuities, Shocks, etc. As I will show In Volume 30 the C-O-R 
constructions of Generahzed Functions enable one to define 'differential 
forms with generahzed coefficients', thus preparing the ground for the 
material in this Chapter serving as foundation for Mechanics with Singular 
Data, the Theory of Splines on nonlinear manifolds, etc. Further, wh en 
combined with the Computational Methods under development at the AI 
Lab of MIT by Gerry Sussman and co-workers this material will be useful 
in the development of Air Traffic Control methodology. 

Another gaol of my work is to develop a general structure for ODE 
systems, to be used in both 'smooth' and 'gene raliz ed' (in the sense of 
Colombeau, Oberguggenberger and Rosinger) Mechanics, Control and 
Numerical Analysis. Since Martin, Crouch have shown that,~in-the linear 
case, Splines may be constructed from linear control system so attention 
will, in the future, focus on the Splines associated with Gene raliz ed Inpouts 
to Nonlinear Control Systems. Work of Sastry and Montgomery indicates 
that imporant examples of such systems will be the Left-Invariant Control 
Systems on Lie Groups, which have been much studied in recent years by 
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Constrained Mechanics and Variational Probl 



researchers interested in Integrable Systems, Robotics, and Aircraft 
Guidance, 



1. Introduction. 

There has been considerable interest recently in constrained 
mechanics and variational problems. This is in part due to Applied 
interests (such as 'non-holonomic mechanics in robotics') and in other part 
due to the fact that several schools of 'pure' mathematics have found that 
this classical subject is of importance for what they are trying to do. I 
have made various attempts [2, 3, 4, 5, 6, 8, 11, 15, 20, 26, 27] at 
developing these subjects since my Lincoln lab days of the late 1950's. In 
this Chapter, I will sketch a Unified point of view, using Cartan's approach 
with differential forms. This has the advantage from the C-O-R viewpoint 
being developed in this Volume that the extension from 'smooth' to 
'generalized' data is very systematic and algebraic. (I will only deal with 
the 'smooth' point of view in this Chapter; I will develop the 'generalized 
function' material at a later point.) The material presented briefly here 
about Variational Calculus and Constrained Mechanics can be found in 
more detail in my boooks, "Differential Geometry and the Calculus of 
Variations" "Lie Algebras and Quantum Mechanics", and "Geometry, Physics 
and Systems". 

Here is the basic set-up. Suppose given the following data: 

A smooth paracompact manifold X (1.1) 

T(X) is its tangent vector bundle (1.2) 

A set {6, co 1 , ..., com] 0 f smooth 1-forms on X. ( 1.3) 

0 is called the action form, {co 1 , ..., co 111 } are the constraint forms. 


Let us suppose given a curve in X: 
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x= {t > x(t) e X: a s t js b}}: [a, b] — > X (1.4) 

dx/dt = v = {t-> dx/dt(t) eT(X):astsb}}: [a, b] ~> T(X) (1.5) 
is its tangent vector or velocity curve. 

(In this Chapter, I suppose all such curves are also smooth.) 

Definition. The following real number associated to the curve 1.3 is called 
the action: 


a(x) =/[a,b]0([dx/dt](t))dt (1.6) 

The following field of 1-covectors along the curve 1.4 is called the force: 


[t~> [dx/dt] (t))J d0} (1.7) 

1.6 and 1.7 are the basic data for both 'mechanics' and 'variational 
calculus'. 

Now, let us deal with 'constraints': 

Definition. The curve 1.4 satisfies the constraints associated with 
the 1-forms {to 1 , ..., co m } iff. it satisfies the following set of Pfaffian 
differential equations: 

0 = to 1 (dx/dt) = ... = co m ( dx/dt) (1.8) 

I will show how the basic Equations of Mechanics can be described 
very compactly and elegantly in terms of this data. 

2* The First Variation formula and the Cauchy characteristic 
curves of de. 
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Keep the data of Section 1. Let us suppose that only the 'action' form 
0 is give, without any constraints. Let x be a smooth curve in X, given as in 
1.4. Suppose that 's', 0 Is s * 1, is a deformation parameter and that 
{s — > xs : [a, b] — > X} is a smooth one-parameter family of curves in X, 
reducing to the given curve x at 's=0'. For t e [a, b], set: 

v(t) = tangent vector to the curve {s --> Xs(t)} at ’t= O' 

The field {t --> v(t) eX X ( t )} of tangent vectors is called an infinitesimal 
deformation of the curve x. Then: 

d(a(Xs))/ ds| s=0 (2.1) 

is called the First Variation of the action function function 1.6 
along the curve x pointing in the direction of the vector field v. 
Using the formula 1.6 for the Action, the Cartan Family Identity: 

'V(e) = VJ de + d(VJ e)': between a differential form and a vector field, and 
Integration-By-Parts, we have the First Variation Formula: 

d(a(xs))/ds| s =o =/[a, b] -[dx/dtj de](v(t))dt + 0(v)(b) - 0(v)(a) (2.2) 


Remark. This formula is a variant of a General Principle: 

The Variational Derivative of the Action is the Force (2.3) 
It also suggests the following: 

Definition. A curve {x: t — > x(t)} is called a Cauchy characteristic 
curve for the 2 -form d0 iff: 

[dx/dt] J d0 = 0. (2.4) 

If x satisfies 2.4, then the First Variation 2.2 vanishes for any infinitesimal 
deformation v which satisifies the following conditions: 
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0(v)(b) = 0(v)(a) = 0 (2.5) 

Conditions 2.5 are called Transversality Conditions. 

At this point, 'symplectic structures and foliations’, 'Hamilton's and 
Lagrange's Equations' (for special choices of 0), etc. enter in a very natural 
way. See [2, 4, 6, 8, 11, 20, 27, 29]. 

3. The differential equations of constr ain ed extrema and the 
augmented action form. 

Let us now suppose that {6, to 1 , .... co m ] is given, as in 1.3. Introduce 
Lagrange Multiplier Variables: 

{Xi, ..., Xxn] (3.1) 

Consider them as Cartesian coordinates of a copy of R m . On X x R m , 
introduce the following augmented 1-form: 

0aug — 0 + Xicul + ... + (3.2) 

Then, 

d0aug = d0 + dX^Aco^ + ... + dXniAo> m + (3.3) 

Xjdtol + ... + Xmdo/ 11 


Definition. A curve {t — > x(t)} in X is an extremal of the constrained 
variational problem associated with the differential form data 1.3 if and 
only if there is a curve in X x R m of the form 

{t -> (x(t), Xi(t), ..., X m (t)} which is a Cauchy characteristic curve of d0 aug . 
In other words, the 'extremas' are the images under the Cartesian 
projection map: {X x R m -> X] of the Cauchy characteristic curves of de au g. 

Theorem 3.1. A curve {t --> (x(t), Xi(t), ..., yt)} is a Cauchy 
characteristic curve for the 2-form 'd0 a ug' if and only if the following 
conditions are satisfied: 


0 = co^dx/dt) = ... = co^dx/dt) 


(3.4) 
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dx/dtj d9 - - Xi(t)[dx/ dt] J dco 1 . ... . X m (t)[dx/dt] J dto 111 (3.5) 

- [dXi(t)/dt] col . ... . [dX m (t)/dt] tom 

Proof. Let v be a tangent vector to the manifold X x R m . Then: 
vj d0 aU g = Vj (de + dXiAtol + ... + dX m Ato m + Xidco 1 + ... + Xindto 131 ) 

* V J d0 + v(Xi)o)l + ... + v(Xi)o)l - (ol(v)dXi - ... - (^(vJdXm (3.6) 

+ Xi[vJ dto 1 ] + ... + X m [vJ do) 01 ] 

3.6 involves one-forms on X x R m . Notice that the only terms on the right 
hand side of 3.6 which involves fdXi, ... dx m } are the terms 
a)!(v)dXi - ... - co m (v)dX m ’. If the tangent vector v is to be Cauchy 
characteristic these forms must vanish. This leads to the condition 3.4. The 

conditions 3.5 now follow from inserting 3.4 into the Cauchy characteristic 
conditions V J d0 aug = O' and using 3.6. 

q.e.d. 

Remark. This result expands the treatment to the 'constrained' case that 
Cartan gave for the 'unconstrained' variational problem in "Lecons sur les 
Invariants Integraux". See [2] for the connection with the traditional 
Lagrange Variational Problem' as expounded in Caratheodory's book and 

for the definition and properties of 'symplectic foliations' and further 
detail. 

The differential equations of constrained mechanics. 

There is considerable confusion in the Literature beween the 
Lagrange Variational Problem (or 'constrained extrema') and 'constrained' 
(and non-holonomic’) mechanics'. I will now describe the latter. Suppose 
again given the following data: 


A smooth paracompact manifold X 


(4.1) 
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T(X) is its tangent vector bundle (4.2) 

A set {e, col, ... ? 0,111} 0 f smooth 1-forms on X. (4.3) 

Definition. Let {x: t --> x(t)} be a curve in X. It is said to be a trajectory 
of the constrained mechanical system associated with the data 
4. 1-4.3 iff. the following conditions are satisfied: 

0 = «)l(dx/dt) = ... = co m (dx/dt) (4.4) 

There is a curve in X x R m of the form 
(t --> (x(t), M-i(t), ..., Mm(t)} such that: 

[dx/dt] J de = m(t)o)i + _ + ^(t) ^ (4.5) 

In other words, 4.4-4. 5 define an ODE system whose solutions are curves in 
X x R m . The 'constrained mechanics trajectories' are the projections in X 
under the Cartesian map projection [X x Rm -> Xj of the solution curves of 
the ODE system 4.4-4.5. 

5. 'Holonomic' constraints. Equivalence of the Constrained 
Extremal and Mechanics equations in the 'holonomic' case. 

Suppose given the following data: 

A smooth paracompact manifold X (5.1) 

A set {0, o)l, ..., «) m } of smooth 1-forms on X. (5.2) 

Indices 1 s a, b, ... s m (5.3) 

Definition. The constraint forms fuA} are said to be holonomic- iff. there 
is a matrix |co a b} of 1-forms such that: 


dco a = o) a b A o) b 


(5.4) 
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Remark. Locally, condition 5.4 is equivalent to the following, more 
'geometric', condition: 

The Pfaffian System {co a = 0} is Frobenius Integrable (5.5) 

Let us now combine conditions 5.4 and the Constrained Extremal 
equations 3.5. The following equations result: 

dx/dtj d0 = - 2 a bMt)[dx/dt] j (cu a bA co b ) - 2 a [dX a (t)/dt] cu a (5.6) 

Rewrite this as follows: 

dx/dtj de = - 2abXa(t)[(a, a b(dx/dt) cub) . oob(dx/dt) cu a b (5.7) 

- 2 a [dXa(t)/ dt] cu a 

The second term on the right hand side of 5.7 vanishes as a consequence of 
the Constraint Equations 4.4, resulting in the following: 

[dx/dt] J de = - 2ab^a(t)[(cu a b(dx/dt) cob) - 2 a [dXa(t)/dt] co a (5.7) 

Theorem 5.1. Let 5.4 be satisfied and let the curve {t --> x(t)} be a 
solution of the Constrained Extremal Equations. Then, {t --> x(t)} is also a 
solution of the Constrained Mechanics Equations 4.4-4.5. 

Proof. That functions ft — > n a (t)} exist satisfying 4.5 is evident from 5.7. 

q.e.d. 


Here is the converse: 

Theorem 5.2. Let 5.4 be satisfied and let the curve {t --> x(t) } be a 
solution of the Constrained Mechanics Equations 4.4-4.5. Then, {t — > x(t)| 
is also a solution of the Constrained Extremal Equations. 
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Proof. We must show that the existence of functions {t --> X a ( t) } satisfying 
5.7 is a consequence of the existence of functions {t --> n a (t)} satisfying 
4.4-4.5. Examining the right hand side of 5.7, we see that the {n a (t) } can be 
obtained by solving an ODE whose coefficients depend on the fx a (t)}. 

q.e.d. 

6. The constrained mechanics equations in a 'Hamiltonian' form. 

So far, we have been working in the context of general manifold 
theory. Let us specialize now to the situation which is close to the 
'Hamiltonian' formalism in the traditional particle mec hani cs case. 


Suppose given the following data: 

n is an integer (6.1) 

The following range of indices: (6.2) 

1 s i, j, ... $n 

X = R 2n +1 = Rn x Rn x R (6.3) 

{q 1 , Pi, t} are Cartesian coordinates on X. (6.4) 

{(q, P, t) --> H(q, p, t)} is a smooth real-valued (6.5) 

function on X, called the Hamilt onian 

0 = 2jpidqi - Hdt (6.6) 

dH = ZiHidqi + 2iH‘dpi + H t dt, (6.7) 


where {Hi, H 1 , H t } are the partial derivatives" of 
the Hamiltonian function with respect to the 
'canonical' coordinates 6.4. 
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Theorem 6.1. de - 2i(dqi - H*dt)A(dpi + Hidt) (6.8) 

Proof. Follows from 6.7 and 6.6, by a direct compilation, which is left to te 
reader. 

Theorem 6.2. Let V be a smooth vector field on X. Then: 

VJ d0 = Zi(V(qi) - H*V(t))(dpi + Hidt) - E t (V(pi) + HiV(t))(dqi - Hkit) (6.9) 
In particular, if: 

V(t) = 1 (6.10) 

then: 

VJ de = Zi(V(qi) - tf) (dpi + Hidt) - £(V(pi) + Hi)(dqi - Hkit) (6.11) 

Proof. Apply the operation 'VJ ' to both sides of 6.8. 6.11 follows from 
substituting 6.10 into 6.9. 

Theorem 6.3. Keep the hyopotheses of Theorem 6.2 and condition 6.10. 
Suppose further that: 

VJ de = ^(tjZiaidq 1 (6.12) 

where {aO are smooth functions on X and {t --> jx( t) } is a real-valued 
function of 't\ Then, the following relations must be satisfied: 

(6.13) 

(6.14) 


V(qi) = Hi 
V( Pi ) + Hi = px(t)ai 

2i(V(pi) + HOH* = 0 


(6.15. 
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Proof. 6.13-.15 results from combining 6.11 and 6.12, and comparing 
coefficients of independent coordinate differentials on both sides of the 
resulting differential form relation. 

Theorem 6.4. Let V be the vector field on X defined by 6.10 and 6.12. 
Then, the orbit curves {t --> (q(t), p(t), t){ of V are solutions of the 
following ODE's: 


dqVdt = dH/dpi 

(6.16) 

dpi/dt = - dH/dqi + n(t)ai 

(6.17) 

2iai(p(t), q(t), t)[dqVdt] =0 

(6.18) 


Proof. Follows from 6.13-6.15. 

Remark. Equations 6.16-6.18 form an ODE system of (2n+l) equations for 
the (2n+l) 'unknowns: {pi( t), qHt), ji(t)}. They are the Hamiltonian 
version of the Lagrange Equations of Motion for Constrained 
Mechanics. (In this case, there is only one 'constraint, namely 6.18. The 
case of more constraints can be handled similiarly.) 

7. Final remarks about generalizations. 

The material in Section 6 suggests a Generalization of material about 
Symplectic Manifolds, Geometric Quantization, etc. from the traditional case 
abstracted from Particle Mechanics (as in the work of Dirac, van Hove, 
Segal, Kostant, Souriau, etc) to a abstract sitaution paralleling the material 
developed in Section 6. 

I will briefly sketch such generalizations. Instead of the 'R^n+i' 
situation of Section 6, suppose that we are on a manifold X, with the 
following relation: 


de = Q - dHAdt 


(7.1) 
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'H' and 't' are smooth functions on X. 0 and Q are, respectively, a 1- and 2- 
form on X and Q is closed. Suppose that Vh is a vector field on X such that: 

VrJ de = (aoj (7.2) 

V H (t) - 1, (7.3) 

where V is a 1-form defining the constraints and V is a function on X. 

7. 1-7.3 imply: 

V H J Q - V(H)dt + dH = nco (7.4) 

This relation generalizes the duality relation between 'infinitesimal 
symmetries' and 'conserved functions' that plays the basic role in the 
'geometric quamtization' theory of unconstrained conservative mechanical 
systems. I plan to study this Geometric Structure further in a later 
Volume. 
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Abstract 


In this report the problem we axe going to study, is the interpolation of a set 
of points in the plane with the use of control theory. We will discover how 
different systems generate different kinds of splines, cubic and exponential, 
and investigate the effect that the different systems have on the tracking 
problem. Actually we will see that the important parameters will be the two 
eigenvalues of the control matrix. 
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1 Introduction 

I would like to begin by thanking my advisor Pr. Anders Lindquist for 
initiating contact with Texas Tech. I would also like to thank Pr. Clyde F. 
Martin for being my advisor at Texas Tech. 

In this report we will look at a way to store signatures. We want to 
do this by storing only a minimal amount of points on the signature curve, 
and still be able to reconstruct the curve by interpolating these points. The 
interpolation will be performed by splines, and we will look at the common 
splines-problem from the control theory point of view. We can construct a 
trajectory of a system that passes through a specified set of points, and thus 
interpolate the points. 

Two questions that need to be answered arise. First, when is it possible 
for the system to pass through the points? Second, when there are many 
ways to accomplish that, what sort of conditions should we demand that the 
system fulfill in order that we get a unique solution. 

The question of when it is possible to interpolate the points will be an- 
swered in the general case in section 2 Reachability and for our particular 
system in section 3 The System . 

An algorithm to find the solution is developed in section 4 Derivations. 

The choice of boundary conditions is discussed in section 5 Boundary 
conditions. 

In section 6 Results the results of tests done upon parametric curves are 
displayed and discussed. 

A summary in Swedish is provided in section 7 Resume - Summarnj in 
Swedish . 

The programs I have been using are included in section 8 Programs . In- 
cluded among the Matlab programs is an altered version of the original Mat- 
lab program quadS. 

When we have answered the two questions, we have to decide what kind 
of system we will use for the interpolation. We can easily imagine that we 
would get to completely different curves if we asked a pedestrian to walk 
through a set of points and if we asked a cyclist to ride his bike through the 
same points. In the first case we would get (if we suppose that man is lazy), 
linear interpolation, and in the second case we would get a smooth rounded 
curve. This is the same as in our case where we have exponential and cubic 
parametric splines. 
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2 Reachability 

In this section we will determine under what circumstances it is possible to 
take a time invariant linear system from a point Xo at time to to a point Xi 
at time It is a vital property to us, because, in order to interpolate we 
have to be able to pass through the points. We will call a system completely 
reachable if it has the property that this can be done in any positive time 
for any two points. 

This is a classical control theory question, and it is answered by the 
following well known theorem, which was at least implicitly discovered by 
Kalman. 


Theorem 2.1 Suppose that the system below is given, 

| x = Ax + Bu 

1 y = Cx 


( 2 . 1 ) 


where A is n x n and B is n x k. Then it is completely reachable iff T = 
[B, AB, A 2 B , . . . , A n ~ l B], is full rank. 

In order to prove this and understand how the reachability concept can be 
characterized by the matrix T, we will have to look at the general solution 
of equation 2.1. 


x{t) = e-^-^Xo + [‘ e A ^~ 3) Bu k (s)ds (2.2) 

JtQ 

In order to have the desired state Xi at time t x the following equality must 
be satisfied. 

Xi = e^-^Xo + r e A{tl - 3) Bu k (3)c[s (2.3) 

JtQ 

The question of reachability is now easily seen to be the question of whether 
there are any solutions to the mapping L :U >->R n such that 

Lu 4 /“ e A{tl ~ 3] Bu(s)ds = Xi - e^-^Xo = d (2.4) 

JtQ 

Since we recognize L as a linear operator, it is as always very fruitful to use 
a theorem from the general theory of functioned analysis [4, p.250]. 
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Theorem 2.2 If X, Y are complete inner-product spaces and A:X » —* Y is 
a linear continuous operator then 

JmA = amAA' 

We know that R n is a Hilbert space, but we have to look at what kind of 
space U is. We choose to introduce the following inner product for U 

(u ,v)u = [ * uWMO* 

Jt 0 

and it can be checked that U becomes a Hilbert space. We know that there 
is a adjoint operator L m : R n >— >IA such that 

(d, Iu)rh = (X‘d,u)w 

and we get the adjoint L m of the mapping L through this equation 

(d,Iu) 8 n = d' / e A(tl ~ 3 ) Bu(s)ds 
J to 

= ^ ^B'e A '^ l ~ 3 u.(s)ds = (L~d,u)u 

We thus have a linear mapping W = LL m : R. n i— *R n , that is, it is actually 
only a matrix operator. 

w = LL' = r e^'-^BB'e-^'-’Us 

Jt 0 

With only the basic knowledge about matrices we will now be able to prove 
the following lemma 

Lemma 2.1 Let A be n x n and B be n X k. Then, for all to, ti such that 
t 0 < t\ we have 

am W{t Q ,t x ) = am [b, AB,A 2 B,...,A n ~ l B] 

Proof: We will do this by showing that amr C 3mW and amW C amr. 

pmr c am w] 
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Let a £ &zvW, which implies that 0 = a'W^a = ff* a'e^ tl ^acis, 

i.e. 



from which it follows that 

flV V(tr,) a5 0, Vae[to,fi], 

i.e, 

E 4f(t, - sYB\A')U = 0. 
i=o I * 

This implies that [i?, AS. A 2 5, . . . , A n ~ x B\ a = 0. That is, for an arbitrary 
a £ RtxW we have a £ £erF' which implies that £er W C £err' and by a 
theorem from fundamental algebra this equals CTmT C Jm W. 

pm W C 3mT] 

Suppose a £ 3mlV. Then there exists a x £ R n such that a = Wx, and 
hence 

a = f; A J B P ^(tx - syB'e A,{tl - 3) xds 

j= o ^ to J- 

from which it is obvious that a € 3m [B, AB, A 2 B , . . .] and by Cayley- 
Haxniltons theorem that a 6 3mT, which concludes the proof. □ 

The main theorem of this section will follow as a direct consequence of 
the lemma. 

Proof:[Theorem 2.1] By lemma 2.1, Xi — e' 4 ^ 1- ‘°^Xo — d € R" and 3mF = X n 
implies complete reachability. □ 
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3 The System 

We will consider the system with the state and dynamics given by 


x = 


x 

X 

y 

L y J 


where 


x = /lx + Bn 
y = Cx 


(3.5) 



° 

1 

0 

0 \ 


° 

0 \ 


0 

Ai 

<*1 

a 2 


1 

0 

A = 

0 

0 

0 

1 

,B = 

0 

0 


,C = 


10 0 0 
0 0 10 



(3.6) 


This gives us the property 

y = Cx = 

and the system dynamics will look like 

f x = \\i + a^y -r a 2 y + ui 

y = A 2 j) + + 02 % + u 2 

And it gives us the following T 

'0010 Aj a 2 

1 0 Aj a 2 &2P2 + Aj Qi + a 2 (A a + A 2 ) r 2>7 r 2 ,s 
r = 

oooi 02 A 2 r 3 , 7 r 3l3 

.0 1 02 A 2 01 + 02^1 + A 2 ) &202 + Aj r 4 ,7 r 4 , 8 


(3.7) 


Tu r lt8 


where 
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Tl.7 



IYs 

= 

OC\ + 0?2(Ax + A 2 ) 

T 2 , 7 

- 

cl\'32 4“ 4“ Q f 2/^2(2Ax A Ao) -f* A^ 

r 2,8 

= 

dVS'i 4“ ai(A x 4“ A 2 ) 4* Aj 4~ ^2 AiA 2 t o 2 A^ 

r 3 , 7 

— 

/3\ 4- /3 2 ( Ax 4- A 2 ) 

r 3 , 3 


a 2/ J 2 4- A; 

r^,7 

= 

o 2 /?2 4* (Ax 4- Ao) 4- # 2 Ax 4- ( -3 2 A; A 2 4" I? A 2 

fVs 


ctx 4“ c*-20\ 4* &2$2 ( A i + 2A 2 ) 4* Ao 


As T is easily seen to have full row rank, by simply looking at the first four 
col umn s. The class of systems we are going to consider in this article will all 
have the desired property of complete reachability, by Theorem 2.1. 


4 Derivations 

Given a set of points in the plane {(z 0 , y 0 ), (xx, yi), . . . , (x n ,yn)} and the 
corresponding time points {to, tu ■ - • , tn} we would like to find the control 
functions {u 0 , Ui, . . . , u„_i } that takes the system through the points at the 
specified times. 

Let’s study the control u* : 

As t 6 [tk,tit+ 1 ] the state of the system will be 

x(t) = e A{t ~ tk) x k + f e A(t ~ 3) Bu k (s)ds (4.3) 

Jt k 

and as we want the state of the system to be x*+x time 1 we get the 
following condition. 

X * + 1 = e' 4(t * +1_f * ) Xfc + / + e /1(f ‘ +1_j) 5ujt(s)ds 

Jtk 



(4.9) 
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The solution u* to equation 4.9 that minimizes the norm of the control 
signal is then given by 

u fe (f) = B'e- A,i e- A3 BB'e- A ''ds^ (e- Atk+i x k+l - e~ Atk x^) (4.10) 

The control would be specified completely by equation 4.10 if we knew the 
whole state- vector at each interpolation point. We know all but we 

do not know the (**,£*)• To determine the 2 (n + 1) unknowns we have to 
apply some conditions on the solution, and our first choice will be to require 
that the control is continuous, 

Assumption 4.1 

u*(£fc+ 1 ) = k = 0 . . . n — 2 


This will give us 2(n — 1) conditions and will leave only four unknowns. We 
will apply the additional conditions on the boundary and we will have to 
come back to this in section 5 Boundary Conditions . 

In many applications it is just the shape of the trajectory that matters, 
and not the velocity that the system tracks the trajectory with. In these 
cases it makes it much easier to assume that the time between each point is 
a constant. 

Assumption 4.2 Let x — i*. = h. 

Assumption 4.2 can be used to simplify the integral in equation 4.10. 

/ U+l e- A3 BB'e- A ' s ds = {r = s - t k } = 

Jtk 

e~ Atk ^ e~ Ar B B'e- A ' r dr e~ A ' tk 
Jo 

V ' 

matrixconstant 

M= [\- Ar BB't- A ' T dT 

Jo 


Definition 4.1 
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We can now rewrite expression 4.10 as 

Uk(t) = B t e- A 'l*- t '>M- 1 {e- Ah x k+ i-Xk) (4.11) 

Using this expression we will now investigate how the continuity condition 
in assumption 4.1 looks. 

u fc (tfc+i) = B / e- 4 'SV/- l (e- AA x fc+l - x*) = 

B , M~ l (e~ Ah Xk+2 — 'Xk+ i) = u fc+ i(t fc+ i) 

We can rewrite this condition using 
Definition 4.2 

Z = M~ l e~ Ah 
W = e~ Ah M~ 1 e~ Ah + M~ l 
giving the equation the simple form 


B' (Z'xk - Wx. k+l + Zxi+ 2 ) = 0, 


In block diagonal form 


Jfc = 0...n — 2 


(4.12) 


■ Z' -w z 
0 Z' -w 


0 0 0 
0 0 0 


Xo 

Xi 

x 2 


0 0 0 
0 0 0 


-W Z 0 

Z' -w z 


= 0 (4.13) 


Xn-2 
Xn-1 
Xn . 


Now, we have to look at what the unknowns are. The vector in equation 4.13 
is made up of subvectors 

‘ fc ’ 

_ * k 
Xk ~ Vk 
. iik . 
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consisting of two knowns and two unknowns. By partitioning the submatrixes 
W, Z as 



and using the notations given in the following the definition, we can keep 
the unkn owns on the left side and move the position coordinates over to the 
right hand side. 

Definition 4.3 
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And we get the system 








r tel - 

*0 

' Zfi 

-w 3 

z? u ... 

0 

0 

0 


v-uei 

X 1 

0 

z?i 

-w 3 ... 

0 

0 

0 


xrvel 

Ao 

0 

0 

0 

-w? 

7 s 

"/ u 

0 


X n-2 

0 

0 

0 

Z,f 

-W 3 

Zt* 


^vel 
X n- 1 






i 

X 

i 








r 1 

‘ z 3 

-w 3 

Zrl ... 

0 

0 

0 


vP° s 

X 1 

0 

Zfi 

-w 3 ... 

0 

0 

0 


X 2 

0 

0 

0 

-w 3 

7 a 

0 


X n— 2 

0 

0 

0 

7 3 

-w 3 

Z B _ 


X n- 1 







Y poj 

L -Si J 


As we evaluate the right hand side, we get a constant vector. Depending 
on what kind of boundary conditions we choose to use, the derivations differ 
from here. We will deal with the most common cases, each case in turn, 
beginning in the next chapter. 


5 Boundary Conditions 

In order to get a unique solution to our problem, we had to apply the con- 
tinuity condition and the boundary conditions. The continuity condition is 
a rather natural condition, but the boundary conditions have to be studied 
more extensively. The four most common choices of boundary conditions in 
the one dimensional case according to [2] are 

1. Zero velocity at the first and at the last point. 

2. Specified starting and ending direction. 

3. Natural boundary conditions, y n = 0. 

4. Periodical conditions. 
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We will look at how the two dimensional equivalent of choice number 1 
and 2 affect the curves, and for which the same derivation is valid. 

Item 3:s two dimensional equivalent, x = y = 0, requires a derivation and 
will be studied in subsection 5.2. 

We will avoid dealing with condition 4 as it only complicates the calcula- 
tions, and would only be natural and interesting if we were to write the same 
word twice, connected with itself as RogerRoger. 

Because the effect of the boundary conditions are similar at both bound- 
aries we will restrict ourselves to only talking about the starting point. 

5.1 Known velocities at boundary 

We have assumed that we also know and xj^. Moving these over to the 
right hand side, we get a block diagonal matrix system to solve. 

: (5.14) 

fin— 2 
ftn-1 

Where 

n, = -W + if’xr-W-W 

fit = -ZjxJT, + - z® X&* 

O _ , U/B v pos _ yB poj _ vB vel 

i‘n— 1 — ^rl^n—2 > ”r X n _ l ^'ru^’n ^fu^n 

This can easily be solved, and with a linear increase of time. Having solved 
the svstem above we now know all the states of the system in the interpolation 
points. We will now use equation 4.11 for the control, and insert it into 
equation 4.8 to get the trajectory. 

x(0 = (x* + t~ Ar B B' e~ A ' T dr - x*)) 

for k = 0 ... n — 1 
t £ [*fc>£fc+i] 


-W, B Zfc ... 0 0 1 [ x“ e ' 

Zfi -WP ... 0 0 xf 

o 0 ... -W? Z B x^ 2 

0 0 ... Zg -w? \ [ x"tx 


(5.15) 
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As we can see from equation 5.15 the fundamental matrix and the integral 
is the same for all k and only has to be evaluated for t — tk between 0 and h, 
once and for all. 

5.1.1 Zero initial velocity 
We can choose to set 
Assumption 5.1 

(i 0 ,y 0 ) = 0 
(in.yn) = 0 


With this condition, we will let the system start up in whatever direction 
that minimizes the energy norm of the control signal and takes the system 
to the second point. 

As we know from one dimensional control theory, a system with a transfer 
function with zeros in the numerator will start off in the opposite direction 
to where it is going. Such undesired properties should certainly be avoided in 
our tracking problem. In the case where c*i = = — $2 = 0? the states 

x and v axe independent, yielding two one dimensional transfer functions. It 
is easily seen to have no zeros, which is good. 

Otherwise, we will have to look at the two dimensional transfer function 
given below: 

o 

_1 

s 4 — (Ai + A 2 )s 3 + (AiA 2 — a 2 /? 2 )s 2 — (ci\pi + a 2 /?i)s — Qt\p\ 


x 


5(5 — A 2 ) cci + a 2 s 
Pi + P2 3 5(5 — Aj) 


( 5 . 16 ) 


This is a bit more tricky, and we will have to find the Q and D of least degree 
that is a solution to the equation 


T(s) = Q(s)D(s)~ l 


(5.17) 


and satisfies 


X(s)Q(s) + Y(s)D(s) = I 


(5.13) 
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It is easy to verify that the choice 


Q(s) = /, D{s) 


s(s — Ax) — (c*i + 
“(j$1 + ,#2s) s(s — A 2 ) 


is a solution to equation 5.17, and there exists X(s) and F(s) so that equa- 
tion 5.18 is satisfied. From Q(s) we can see that the system has no zeros. 
The zero initial conditions should thus not give us any problems. 


5.1.2 Derivative approximation 

Instead of setting the velocity equal to zero, another alternative is of course to 
specify a starting velocity. However, this requires that we make a good choice, 
to avoid situations as exemplified below. Using a bad direction and a high 
velocity boundary condition on y = x~, we get the graph of figure 1. Even 
as we are using n=40 to reproduce the curve, the bad boundary conditions 
axe still ruining the tracking. 



Figure 1: Trajectory of system tracking y — x~, with badly specified starting 
direction 

As discussed in [3, p.86], the fact is that when we set the boundary 
conditions in the parametric case, we do not only specify a direction, but 
also the speed in that direction. The greater the speed, the greater impact 
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the boundary conditions will have on the solution. We axe thus forced to 
make a good choice, and we would like the choice to get better the more 
interpolation points we are using. A simple choice that satisfies this is 

Assumption 5.2 

(i o,yo) 

(in, y n ) 


X! - Xo 

h 

X n X n _i 

h 


which imply that we will set off from the first point in the direction of the 
second point, and arrive at the last point in the direction from the next last 
point. 

Another way of deciding the initial velocity would be to use the same 
technic as in Bezier curves and choose the settings graphically. This would 
probably be the best way to get the desired properties of the signatures. As 
discussed in [3] the choice of boundary conditions will affect the whole curve, 
and the solving of the blockdiagonal system must be done over from scratch, 
making this method a bit slow. If we are going to do this only once to store 
a signature it does not matter. What matters with this method is that it 
adds four more parameters to be stored, and we could probably get equally 
good results just by increasing the number of interpolation points by two. 


5.2 Constant velocity 

Suppose we want to use the boundary conditions 
Assumption 5.3 

(i 0 ,z/o) = 0 

(*^n > Vn ) 0 


This will let the initial direction and constant velocity of the system be 
decided so that the control energy is minimized. Using the system dynamics 
equations 3.7, we get the equation system below. 
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( ft 2 + w < B - «») <“ - z «- = ( 3 - 23) 

= (u„ - wf - [ £ ] ) xr + 

Adding these two equations to equation 4.14 yields a blockdiagonal system 
to solve. This system is two blocks bigger than the one in section 5.1, but it 
can also be solved with a linear increase of time. And once it is solved, we 
can still use equation 5.15. 

A comparison between the three boundary conditions, BC=1 zero initial 
velocity, BC=2 derivative approximation, and BC=3 zero initial acceleration 
is made in section 6. 
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6 Results 

The first tests with the algorithm were done with matrices A with both 
eigenvalues equal to zero, Aj = 0, A 2 = 0, and a\ = 012 = ,3\ = 3? — 0. 
This produces cubic parametric splines, and makes the calculation of the 
fundamental matrix easy. The cubic splines produce smooth curves, but were 
also able to reconstruct a cusp much better than you would have guessed at 
first, as shown in figure 2. 



Figure 2: Reconstruction of y = with cubic parametric splines where 
n=10. 


If we look at the function y = we know that this function describes a 
cusp. But if we parameterized it like x = t 3 ,y = t 2 we see that both x and 
y are smooth cubic functions of t, so it is not very remarkable that it can be 
reconstructed well. 

When we used A-matrices with nonzero eigenvalues, and decoupled x 
and y coordinates, i.e. a x = a 2 = Pi = = 0, we were able to generate 

exponential parametric splines with the basis functions 1, Ajf, e"' 1 *, e 1 and 
l,A 2 /,c Aj ‘,e~ A3 ‘ for the x and y coordinates. 

The result of taking big eigenvalues is almost linear interpolation, which 
can be good for certain applications, but not if it is to be used for storing 


signatures. It is quite obvious that the roundness of a persons signature is 
one important factor of it’s characteristics. Therefore, it’s vital that one of 
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data stored on the signature is the eigenvalues of the A- matrix, as can be 
shown in figure 3. 

Figure 4 shows the corresponding control signals. In the At = 1, Aj == 1 
case the linear part of the control is dominating, and in the Ax = 100, A^ = 
100 case the exponential part is dominating. It is also evident that the 
magnitude of the control signals is greater in the case of larger eigenvalues, 
but that is not a problem for our fictitious system. 

In the case of nonzero we get a coupling between the x 

and y coordinates. This will give us very complicated basis functions, like 
polynomials times exponentials times sine- and cosine-functions. 

As with all approximation methods one should always investigate how big 
the errors are. To do this we had to somehow determine the distance between 
the original curve and the interpolating one. The tests were performed on 
known parametric curves, so we had an explicit expression for the points on 
the original curve. We had to try to find the nearest point on the original 
curve to the point on the trajectory. This was solved numerically with the 
“Golden Intersection algorithm' 1 . For the method to work we have to assume 
two things : 

1. The section of the original curve between the two interpolation points 
nearest in time is convex. 

2. The point on the original curve nearest the point on the interpolating 
curve lies on the section in item L 

As an error estimate I have calculated the distance between a number of 
points on the reproduced curve and the original curve and divided with the 
number of points. We have applied this error estimate method on four dif- 
ferent curves and with different number of interpolation points and 40 points 
between each of these. 


points 

BC = 1 

BC = 2 

BC = 3 

n = 10 

6671.1 

5821.1 

4515.9 

n = 20 

1028.4 

765.6 

502.1 

n = 100 

8.3 

5.9 

3.5 


Table 1: \i units of error for unit circle. 
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Figure 5: Graph, of circle, reproduced with n=10 and A a — 10, A 2 — 10, 
BC=1,2,3 


6.1 Four test curves 

The first curve we tested was the unit circle. This very round curve was 
tracked best by the cubic splines, but the exponential splines did a good job 
too, as can be seen in figure 5 and table 1. We can also see that the error 
was smallest in the case of zero initial acceleration boundary conditions. 



Figure 6: Graph of y = |x|, reproduced with n=10, n=20 and A t — 10, A 2 — 
10, BC=2 


Next, we looked at a curve with the opposite properties, linear and non- 
differentiabte, y = |x|. The error is mainly located between the two points 
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points 

BC = 1 

BC 7 = 2 

BC = 3 

n = 10 

3450.3 

3446.3 

3450.3 

n = 20 

972.6 

972.6 

972.5 

n = 100 

40.7 

40.7 

40.7 


Table 2: ft units of error for y = |x|. 


next to the non-differentiable point. As could be guessed the tracking of 
the curve y = |xj got better the bigger the eigenvalues of the A matrix 
were, the error was reduced to 142 ft units for \\ = A? = 500 in the case of 
BC = 1 and n - 10. but for bigger values the numerical calculations failed. 
Using negative eigenvalues gives the same error as the positive, which can be 
expected since we have a symmetric curve and time interval and by looking 
at the basis functions. The results with different boundary conditions were 
very much alike, as seen in table 2. 

This was the only case were BC 3 did not give us the smallest error. 

Can we get a smaller error with any choice of the coupling parameters. 
Yes, for example by choosing a\ = —0C2 = Pi = “^2 = 10, we get the 
error 3249 fi units for n=10. Choosing these parameters could thus be a way 
to reduce the error, but by using n = 12 instead, we got the error 2506 fi 
units. So we do not get a more efficient way to store it, unless we can find 
parameters so we get below 2506 ft units. 


points 

BC = 1 

BC = 2 

BC = 3 

n = 10 

6506.6 

5313.7 

3842.9 

n = 20 

1014.4 

725.4 

460.9 

n = 100 

8.3 

5.3 

3.4 

Table 3: 

/j units of error for cycloid. 

points 

BC = 1 

BC = 2 

BC = 3 

n = 10 

13897.6 

11664.3 

8867.6 

n = 20 

2100.5 

1531.0 

1001.7 

n = 100 

17.7 

11.3 

7.0 


Table 4: ft units of error for prolate cycloid. 







Figure 8: Graph of p 
and X\ = 10, A 2 = 10 
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Finally we look at a cycloid. 

f x = Tt — sin rr t 
1 y = 1 — cos rrt 

and a prolate cycloid. 

f x = ~t — 2 sin xt 
| y - 1 - 2 cos xi 

It is evident that the cusp and the crossover do not cause any problem, as 
could be expected since we are using a parameterized interpolant. 

We can compare the difference of the graphs in figure i and figure 8, 
where we are using one crude approximation with n=10 and one extensive 
aDproximation with n=l00. This time it is evident that the error is m ain ly 
located at the boundaries. In table 3 and table 4 we can see that the best 
results came from using constant velocity boundary conditions. 

Looking at table 3 and table 4 again, we see that the error decreases at 
an approximately cubic rate as the number of interpolation points increase, 
which is much better than the quadratic decrease that can be seen in table 2. 

\Iv guess is that this behavior comes from the fact that the curve y — \x\ 
does not have a differentiable parameterization. 

6.2 Applied on a signature 

Included as figure 9 is a scanned picture of my own signature. I have tried 
to pick some roughly equidistant points on the signature, (According to the 
scale indicated on the axis.) and used the interpolation algorithm we have 
been studying to reproduce it. The reproduction is made with n = / 4, 
x l = A 2 = 10 and no coupling between the two coordinates. For boundary 
conditions I have chosen to use constant velocity, since it has been the most 
successful condition. 

As can be seen we get a very close resemblance between the original and 
the reproduction. How close is hard to say because we do not have the 
signature given as a parameterized function, therefore we are not able to 
calculate the error as before. 

The things characterizing the signatures, are also the things that are hard 
to recover with the interpolation. Such as the turnover in the connection from 
the "P” , and the cusps in the “r : ’ . To get a good reproduction, an equidistant 
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distribution of the interpolation points is not enough, more points has to be 
concentrated around the characterizing areas. 

7 Resume - Summary in Swedish 

Lagring av signaturer kan goras pa manga satt. Vi har valt att lagra ett antai 
punkter pa signaturen, och sedan reproducera denna genom att mterpolera 

Dunkterna med splines. > , 

Genom att anvanda valkanda resultat inom systemteon sa kan man gener- 

era olika sorters splines genom att andra pa nagra parametrar. Jag har nytt- 
jat denna metod for att generera parametriserade splines i planet. Man rnser 
snart att man maste infora randvillkor pa losningen. och valet av dessa far 
inte ske hur som heist eftersom de paverkar hela losningen. Darfor har jag 
lagt ner en hel del arbete just pa denna punkt. De basta resulteten har jag 
erhallit genom valet att ha konstant hastighet vid andpunktema. 

En av de saker som karaktariserar handstilar ar dess rundhet. Denna kan 
ges en direkt Sversattning i egenvardena till systemmatrisen, och vi kommer 
Illtsa valja att lagra dessa utover punktema pa signaturen. 
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8 Programs 

The progr ams in Matlab and Maple that were used to implement the algo- 
rithm developed in this report follow. 

8.1 Matlab Program 

To make it easier to understand the structure of the program, the following 
flow charts describe how the programs are traversed. 



The loops marked with an unfi lled circle is only available when the inter- 
polated points are given by the parameterized function (XX(t),YY(t)). 
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function error=run(BC,nn) ; 

7.BC '/. Type of boundary conditions 

7. 1 = zero initial and final velocity. 

7. 2 = heading for first point, last. 

7. 3 = zero initial and final acceleration. 

7, If nn is specified, runs alg with nn points given by XX, YY. 
7. Otherwise runs alg with points given by ginput. 
global n h alphal aipha2 betal beta2 lambdal lambda2 
global errorcalc ctrlsignal 
clg 


7.7. Setting of parameters 7.7. 

alphal=0 ; 

alpha2=0 ; 

betal=0 ; 

beta2=0 ; 

lambdal=10; 

lambda2=10 ; 

7. Decides what steps are going to be made 
errorcalc=l; 7. error estimation 

ctrlsignal=0; 7. plotting of control signals 


if nargin == 1 
[x , y] =ginput ; 

R ( 1 , : ) =x * ; R(2,:)=y’; 
else 

R=points(nn) ; 
end; 

n=length(R)-l; 7. Number of interpolationpoints -1. 

h=2/n; 7. Time inbetween points 

7.7. Plots a circle at all the points thats interpolated 7.7. 
hold on 
for i=l:n+l 
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plot(R(l,i) ,R(2,i) , ’o’) 
end; 

7.7. Calling alg with, the prepared data 7.7. 
if BC == 3 

error=alg2(R) ; 
else 

if BC == 2, xvelO=(R( : ,2)-R( : , 1) )/h; 

else xvelO=zeros (2 , 1) ; 
end; 

if BC == 2, xveln=(R( : ,n+l)-R( : ,n) )/h; 

else xveln=zeros (2 , 1) ; 
end; 

error=alg(R,xvelO , xveln) ; 
end; 

7.7. Loop to allow graphic alteration of BC. 7.7. 

b=input(’"l" for graphic mod of BC, "0" to quit ’); 
while b == 1 , 

xvel0=10*(ginput (1) ’ -R( : , 1) ) ; 
xveln=-10*(ginput(l) ’~R( : ,n+l) ) ; 

cl S 

hold on 
for i=l:n+l 

plot(R(l ,i) ,R(2,i) , ’o’) 
end; 

error=alg(R,xvelO, xveln) ; 

b=input(’"l" for graphic mod of BC, "0" to quit ’); 
end; 

end; 


function R=points (nn) ; 

7, Forms R with the help of XX, YY. 
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7 , R = 2* Cnn+1) -matrix, 
global a b h 


a=l ; b=a; 
h=2/nn; 

for i= 0:nn 

R( : ,i+l)=[XX(-l+i*h) ; YY(-l+i*h)] ; 
end 
end; 


function res=XX(t) ; 
global a b 

res= a*t*pi-b*sin(t*pi) ; 
end; 


function res=YY(t) 
global a b 

res= a-b*cos(t*pi) ; 
end; 


function error=alg(R,xvelO .xveln) ; 

7. R = Matrix of interpolationpoints , xO ... xn. 

7. first row = x-coordinates , second row = y-coordinates . 

7, xvelO, xveln = Boundary conditions 

global a b A B C n h. alphal alpha2 betal beta2 lambdal lambda2 
global errorcalc ctrlsignal 

7.7. The System 7 , 7 , 

A=[[ 0100] 

[ 0 lambdal alphal alpha2] 
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[ 0 0 0 1 ] 

[ betal beta2 0 lambda2]] ; 

B= [ [0 0] 

Cl 0] 

Co o] 

C0 1]]; 
c=CCi o o o] 

Co 0 1 0]] ; 

7.*/. Calculation of the integral from 0 to h in m steps 7.7. 
m =40; 7. number of points between interpolationpoints 

tol=le-08; 7. the numeric error tolerance 

MtauC : ,l:4)=zeros(4) ; 
tau=0 ; 
for j=l:m 
oldtau=tau; 
tau=oldtau+h/m ; 

Mtau(: ,4*j+l:4*j+4)= quad8mod( ' int ’ ,oldtau,tau,tol) 

+ MtauC : ,4*j-3:4*j) ; 

end; 

M=Mtau(: ,4*m+l :4*m+4) ; 


7.7, Forming of the Matrixes for the Blockdiagonal system 7.7. 


e_Ah=expm(-A*h) ; 
Minv=inv(M) ; 
ZZ=Minv*e_Ah; 
W=e_Ah’*ZZ+Minv; 


WL=Cwy(2,2) WW(2,4); WW(4,2) WW(4,4)] ; 7. Partitioning matrixes 

ZLU=[ZZ(2,2) ZZ(2,4); ZZ(4,2) ZZ(4,4)] ; 

ZLL=[ZZ(2,2) ZZC4.2); ZZC2.4) ZZ(4,4)] ; 

WR=CWW(2,1) WWC2.3); WW(4,1) WW(4,3)]; 

ZRU=[ZZ(2, 1) ZZC2.3); ZZ(4,1) ZZC4.3)] ; 

ZRL=[ZZ(1,2) ZZ(3,2); ZZCl.4) ZZ(3,4)] ; 
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The boundary conditions 7.7. 

xvelC : , l)=rvelO; 
xvel( : ,n+l)=xveln; 

7,7. Forming of the right side of the Blockdiagonal system /./. 
for i=2:n 

Omega( : , i)=ZRL*R( : , i-1) -WR*R( : , i)+ZRU*R( : , i+1) ; 
end; 

Omega ( : ,2)=0mega( : ,2)+ZLL*xvel( : , 1) ; 

Omega ( : ,n)=0mega(: ,n)+ZLU*xvel( : ,n+l); 

7,7. Gausselimination to produce upper triangular system 7.7. 

DD(: , 3:4)=WL; 
for i=3:n 

zd=ZLL*inv(DD( : ,2*i-3:2*i-2)) ; 

DD( : ,2*i-l:2*i)=WL-zd*ZLU; 

Omega ( : ,i)=0mega( : , i)+zd*0mega( : ,i-l) ; 
end 


7.7, Backsubstitution to solve for the xvel 7.7. 

xvel ( : , n) = DD ( : , 2*n- 1 : 2*n) \0mega( : ,n) , 
for i=n-l: ~1 : 2 

xvelC:, i)-DD(:,2*i-l:2*i)\(ZHJ*xvel(:,i+l)+0aega(:,i)); 

end; 

7.7. Making of state vectors 7.7. 
for i=0:n 

x( : , i+l)=[[R(l,i+l)] 

[xvel(l , i+1)] 

[R(2,i+1)] 

[xvel (2, i+1)]] ; 
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function error=alg2(R) ; 

7. R = Matrix of interpolationpoints , xO . . . xn. 

7, first row = x-coordinates , second row = y- coordinates . 
7 , BC = acceleration in x and y direction are both = 0. 


global a b A B C n h alphal alpha2 betal beta2 lambdal lambda2 
global errorcalc ctrlsignal 

7 . 7 . The System 7 . 7 . 

A= [ [ 0 1 0 0] 

[ 0 lambdal alphal alpha2] 

[ 0 0 0 1 ] 

[ betal beta2 0 lambda2]]; 

3= [CO 0] 

Cl o] 

Co o] 

C0 1]]; 

C=[[l 0 0 0] 

[0 0 1 0 ]] ; 

7 . 7 , Calculation of the integral from 0 to h in m steps 7 . 7 . 

m =40; 7. number of points between interpolationpoints 

tol=le-08; 7. the numeric error tolerance 

MtauC : , 1 :4)=zeros(4) ; 
tau=0; 
for j=l:m 
oldtau=tau; 
tau=oldtau+h/m; 

MtauC : , 4*j + l : 4*j+4) = quad8mod( ’ int ’ ,oldtau,tau,tol) 

+ MtauC: ,4*j-3:4*j); 

end ; 

M=MtauC : ,4*m+l : 4*m+4) ; 
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7.7. Forming of the Matrixes for the Blockdiagonal system 7.7. 

e_Ah=expm(-A*h) ; 

Minv=inv(M) ; 

ZZ=Minv»e_Ah; 

W=e_Ah.’*ZZ+Minv; 

WL= [WW (2,2) WW(2 ,4) ; WW(4,2) WW(4,4)]; 7. Partitioning matrixes 

ZLU= [ZZ(2 ,2) ZZ(2,4) ; ZZC4.2) ZZC4.4)] ; 

ZLL=[ZZ(2,2) ZZC4.2); ZZ(2,4) ZZ(4,4)] ; 

WR= [WW (2,1) MWC2.3); WWC4.1) WWC4.3)] ; 

ZP.U= [ZZ (2 , 1 ) ZZC2.3); ZZ(4,1) ZZC4.3)] ; 

ZRL-[ZZ(1,2) ZZ ( 3 , 2) ; ZZ(1,4) ZZ(3,4)] ; 

Ulu= [Minv(2 , 2) Minv(2,4) ;Minv(4,2) Minv(4,4)]; 

Uru= [Minv(2 , 1) Minv(4, 1) ;Minv(2,3) Minv(4,3)]; 

VT=[Iambdal alpha2; beta2 lambda2] ; 

Vr=[0 alphal ; betal 0]; 

7.7, Forming of the right side of the Blockdiagonal system 7.7. 
for i=2:n 

OmegaC: ,i)=ZRL*R( : ,i-l)-WR*R( : , i)+ZRU*R( : ,i+l) ; 
end; 

OmegaC: ,l)»(Vr-Uru)*R(: ,1) + ZRU*RC:,2); 

OmegaC : ,n+l)=ZRL*R( : ,n) - CWR-Uru+Vr)*R( : ,n+l) ; 

77. Gauss elimination to produce upper triangular system 7.7. 

DDC : , 1 : 2)=Ulu-Vl ; 
for i=2:n+l 

zd=ZLL*invCDD(: ,2*i-3 : 2*i-2) ) ; 

DDC : ,2*i-l:2*i)=WL-zd*ZLU; 

OmegaC : ,i)=0megaC : ,i)+zd*0megaC : ,i-l) ; 
end 

DDC : ,2*n+l:2*n+2)= (WL-Ulu+Vl) - zd*ZLU; 
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7.7. Backsubstitution to solve for the rvel 7.7. 

xvel(: ,n+l)=DD( : , 2*n+l : 2*n+2) \0mega( : ,n+l) ; 
for i=n:-l:l 

xvel( : , i) =DD( : , 2*i-l : 2*i) \(ZLU*xvel( : , i+D+CmegaC : ,i)) ; 
end; 

7 . 7 . Making of state vectors 7.7. 
for i=0:n 

x(: ,i+l)-[[R(l,i+l)] 

[xvelCl, i+1)] 

[R(2,i+1)] 

[xvel(2,i+l)3] ; 

end; 


7.7. Plotting of trajectory 

7.7. and error estimate calculation 

sumnorm=0 ; 
hold on 
for j=0:m 

eAtau=expm(A*j*h./m) ; 
for i=0:n-l 

entry=eAtau* (x ( : ,i+l)+Mtau( : ,4*j+l :4*j+4)*Minv* 

(e_Ah*x( : ,i+2)-x( : ,i+l))) ; 
plot (entry (1) ,entry(3) , ’ . J ) 
if errorcalc 

sumnorm = sumnorm + dist(entry(l) ,entry(3) ,i,h) ; 
end; 
end; 
end; 

7.7, Plotting of the control signals 7.7. 
if ctrlsignal 
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pause, clg 

subplot(2, 1 , 1) .hold on 
subplot (2 , 1 , 2) .hold on 
for i=0:n-l 
for j=0:m 

csignvec ( : , j + 1) =3 ’ *expm(-A ' * j *h/ 40) *Minv» 
(e_Ah*x ( : , i+2)-x( : , i+1)) ; 

end; 

subplot (2, 1,1) ,plot(i*h:h/m: (i+l)*h,csignvec(l, :)) 
subplot (2, 1,2) .plot (i*h : h/m: (i*l)*h, csignvec (2, :)) 
end; 
end; 

error=sumnorm/ n/m ; 
end; 


function [Q,cnt] = quad8mod(?,a,b,tol) 

'/.Alteration of the original matlab toolbox program. QUAD 8 
7, Numerical evaluation of an integral, higher order method. 

7. Q = QUAD8(’F’ ,A,B,T0L) approximates the integral of F(X) 

7. from A to B to within a relative error of TOL. 

7. ’ F ’ is a string containing the name of the function. 

7. The function must return a 4*4-matrix output value if 
7. given an input value . 

7. Q = Inf is returned if an excessive recursion level is 
7. reached, indicating a possibly singular integral. 

7, QUAD 8 uses an adaptive recursive Newton Cotes 8 panel rule. 

7. Cleve Moler, 5-08-88. 

7. Copyright (c) 1984-94 by The MathWorks , Inc. 

7. [Q,cnt] = quad8(F,a,b,tol) also returns a function 
7, evaluation count. 

7. Top level initialization, Newton-Cotes weights 

w= [3956 23552 -3712 41984 -18160 41984 -3712 23552 3956] /14175 ; 

x=a + (0:8)*(b-a)/8; 
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7. set up function call 
for i=x 

y = [y fevalCF, i)] ; 
end; 

'!, Adaptive, recursive Newton-Cotes 8 panel quadrature 
QO = zeros (4) ; 

[Q.cnt] = quad8stpmod(F , a,b , tol , 0 , w , x,y , CJ0) ; 

cnt = cnt + 9; 

end; 


function [Q,cnt] = quad8stpmod(F , a,b , tol ,lev, w,xO ,f 0,Q0) 
^Alteration of the original mat lab toolbox program. QUAD8ST? 

7. Recursive function used by QUAD8. 

'!• [Q , cnt] = quad8stp (F , a , b , t ol , lev , w , f , QO) tries to approximate 
7. the integral of f(x) from a to b to within a relative error 
7. of tol. 

7. F is a string containing the name of f . The remaining 
7. arguments are generated by quad8mod or by the recursion. 

% lev is the recursion level. 

7. v is the weights in the 8 panel Newton Cotes formula. 

7, xO is a vector of 9 equally spaced abscissa is the interval. 

7. fO is a matrix of the 9 function values at x. 

7. QO is an approximate value of the integral. 

7. Cleve Moler , 5-08-88. 

7. Copyright (c) 1984-94 by The MathWorks, Inc. 


LZVMAX = 10; 

7. Evaluate function at midpoints 
7. of left and right half intervals, 
x = zeros (1 , 17) ; 
x(l : 2 : 17) = xO; 

x C 2 : 2 : 16) = (x0(l:8) + x0(2:9))/2; 

f C : ,1:4)= fOC : ,1:4); 
for i=l:8 
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f(: ,8*i-3:8*i) = f evalCF ,x(2*i) ) ; 
f ( : ,8*i+l:8*i+4) = fO(: ,4*i+l:4*i+4) ; 
end; 

7, Integrate over half intervals, 
h = (b-a)/16; 

Q1=0 ; Q2=0 ; 

for i=l:9 

Q1 = qi + h*v(i)*f ( : ,4*i-3:4*i) ; 

Q2 = Q2 + h*w(10-i)*f ( : , 69-i*4:72-i*4) ; 
end; 

Q = Ql + Q2; 

7, Recursively refine approximations, 
if normCQ - QO) > tol*norm(Q) k lev <= LEVMAX 
c = (a+b)/2; 

[Ql.cntl] = 

quad8stpmod(F , a, c , to 1/2 , lev+1 , w,x( 1 : 9) ,f ( : ,1:36) , Ql) ; 
[Q2,cnt2] = 

quad8stpmod(F , c , b , to 1/2 , lev+1 , v, x(9 : 17) ,f C : ,33:68) , Q2) ; 
Q = Ql + Q2; 

cnt = cnt + cntl + cnt2; 

end 
end ; 


function res = integrand(v) 
global ABC 

e_AvB=expm(-A*v) *B ; 
res = e_AvB*e_AvB ’ ; 
end; 


function [d] =dist(xx,yy , i ,b) ; 

7, Initiating search algorithm. 
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8.2 Maple Program 

withdinalg) ; 
with (student) ; 

alpha! : =1 ; 
alpha2 : =1 ; 
betal : =0 ; 
beta2 : =0 ; 
lambdal : =100 ; 
lambda2 : =100 ; 

a: =1 ; 
b:=a; 
h: =0 . 2 ; 
n: =22 ; 
m:=15 ; 

R:=vector(n+l) ; . 

XX:=t->a*t*Pi-b*sin(t*Pi) ; 
YY:=t->a-b*cos(t*Pi) ; 
for i from 0 to n 
do 

R[i+1] :=matrix([[XX(-l+i*h)] , 

[YY(-l+i*h)]]) ; 

od; 

A:=matrix( [[ 0, 1, 0, 0], 

[ 0 , lambdal , alphal , alpha2] , 

[ 0, 0, 0, 1] , 

[ betal, beta2, 0, lambda2] ] ) ; 
B : =matrix ([[0,0] , [1,0] , [0,0] , [0,1]]) , 

C : =matrix( [[1,0, 0,0] ,[0,0, 1,0]]), 

alias(Id = &*()) 

Aprim : transpose (A) ; 

Bprim:=transpose(B) ; 

e_At : =t->exponential(-A*t) ; 
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e_Aprimt : =t->erponential(-Aprim*t) ; 
eAt : =t->exponential (A*t ) ; 
e_Ah: =e_At(h) ; 


integranden: 3 proc (v) 

svalm(e_At(v) B k* Bprim k* e_Apnnt(v)) , 

end ; 

integrera: 3 proc (funk) 
global v; 
int (funk, v) ; 
end; 

map (integrera, integranden(v) ) ; 
integralen: =map(simplify , ") ; 

evaluera:=proc (funk) 
global tau; 
subs ( v=tau , funk) ; 
end; 

Mtau:=vector(m+l) ; 

Mtau [1] : =matrix( [ [0 ,0 ,0 ,0] , [0 ,0 , 0 ,0] , [0 ,0 ,0,0] , [0,0 ,0 ,0] ] ) , 
tau: =0 ; 

MO : =evalm(map(evaluera , integralen) ) ; 

for j from 1 to m 
do 

tau:=j*h/m; 

Mtau[j + l] : =e valm (map (evaluera, integralen) -M0) ; 
od; 

M:=Mtau[m+l] ; 

Minv : =evalm(inverse (M) ) ; 

ZZ : =evalm(Minv&*e_Ah.) ; 

WW : =evalm(transpose (e_Ah) &*ZZ+Minv) ; 
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WL: =submatrix(WW, [2,4] , [2,4] ) ; 
ZLU:=submatrix(ZZ, [2,4] , [2,4] ) ; 
ZLL:=transpose(submatrix(ZZ, [2,4] , [2,4] )) ; 

WR : =submatrix(WW , [2 , 4] , [1 , 3] ) ; 

ZRU : =submatrix(ZZ , [2,4] ,[1,3]); 

ZRL: transpose (submatrixCZZ, [1,3] , [2,4] )) ; 


xvel : =vector(n+l) ; 

xvel[l] :=evalm((R[2]-R[l])/b) ; 

xv el [n+l] : =evalm( (R[n+1] -R[n] ) /b) ; 

Omega: =vector (n+1) ; 
for i from 2 to n 
do 

Omega [i] : =ZRL£*R[i-l] -WR&*R[i] +ZRU4*R[i+l] ; 
od; 

Omega [2] :=0mega[2]+ZLL4*xvel[l] ; 

Omega [n] : =0mega[n] +ZLU&*xvel [n+l] ; 

DD : =vector(n+l) ; 

DD [2] : =WL; 

for i from 3 to n 

do 

zd: =evalm(ZLL&*inverse(DD [i-1] ) ) ; 

DD [i] :=WL-zd&*ZLU; 

Omega[i] := 0 mega[i]+zd&* 0 mega[i-l] ; 
od; 

xv el [n] : =linsolve(DD [n] ,0mega[n]); 
for i from n-1 by — 1 to 2 
do 

xvel [i] :=linsolve(DD[i] ,ZLU&*xvel [i+1] +0mega[i] ) ; 
od; 

x:=vector(n+l) ; 
for i from 0 to n 
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do 

x[i+l] : “matrix ( [[R [i+1] [1,1]] , 

[xvel [i+1] [1,1]] , 

[R[i+1] [2,1]] , 

[xvel [i+1] [2,1]]]); 
od; 

plotlist : = [] ; 
for j from 1 to m 
do 

tau:=j*h/m; 

eAtau:=evalm(eAt(tau)) ; 

for i from 0 to n-1 
do 

entry : =evalm(eAtau&*(x [i+1] +Mtau [j+l]4*Minv 

4* ( e_AM*x [i+2] -x [i+1] ) ) ) ; 
plotlist :={ [entry [1 , 1] .entry [3,1]] .op(plotlist)} ; 
od; 
od; 

plot (plotlist, style=point) ; 
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Abstract 


When trying to fly an aircraft as smoothly as possible it is a good idea to 
use the derivatives of the pilot command instead of using the actual control. 
This idea was implemented with splines and control theory, in a system that 
tries to model an aircraft. Computer calculations in Matlab shows that it is 
impossible to receive enough smooth control signals by this way. This is due 
to the fact that the splines not only try to approximate the test function, 
but also its derivatives. A perfect traction is received but we have to pay in 
very peaky control signals and accelerations. 
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2 Introduction 

In the beginning our intention was to calculate control laws for an aircraft 
model so it would fly as smooth as possible in three dimensions. The comfort 
for the passengers was the most important consideration when we forced our 
system, the representation of the plane, to follow a certain trajectory. 

All the programming has been realized in a numeric computation and 
visualization software called Matlab. Also Maple has been used in some 
of the heaviest calculations and my third contact with the more advanced 
computer world was Latex that this report is written in. The second half of 
this paper consists of Matlab code ended with listed references. 

The test began in one dimension with three different kinds of systems. 
By way of introduction the essential conceptions of reachability and stability 
were examined and written down in chapter 3 and 4 respectively. With these 
tools we could investigate the main features of the systems and obtained the 
result that all of them were completely reachable, stable but not guaranteed 
input-output stable (see chapter 5, The Systems). 

A spline is the curve of an n-degree polynomial that is joined in its end- 
points with similar polynomials. They are connected in the way that they 
have the first n-1 derivatives, at the jointly point, in common. Chapter 6 
consists of calculations for the spline approximation and the control theory. 

Chapter 8, Results, discusses some of the results we received and also 
displays examples of graphs that were obtained. The test could not be con- 
cluded in the way we thought due to a surprising combination between con- 
trol theory and splines. The last two pages in chapter 8 deals with this main 
result. 
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3 Reachability 


When controlling an aircraft we will be sure that a suitable control signal u 
can take us to all desirable states. Transfered to our one dimensional case 
we have to determine under what circumstances there is an input signal u 
which transfers the state from x(t 0 ) = Xq to x(tj ) = x t . This is a basic issue 
in systems theory and it leads to the concept of reachability. We will call a 
system completely reachable if it has the property that this can be done in 
any positive time for any two points. 

Most of the trains of thought in the following proofs are derived from 
lecture notes given by Tomas Bjork, Optimization and Systems Theory, Royal 
Institute of Technology, Stockholm, Sweden, during the fall of 93. 

Consider the system. 

x(t) = A(t)x(t) + B(t)u(t ); x(t 0 ) = x 0 - 


with the general solution 

x(t) = $(t,t 0 )x 0 + [ <P(t,s)B(s)u(s)ds. 

Jto 

In order to reach the desired state x(ti) = x\ the following equality must be 
fulfilled. 

rti 

x t = $(t , , t 0 )x 0 + / <P(t 1 ,s)B(s)u(s)ds. 

Jto 

Define d = xj — <P(t t , t 0 )x 0 and let U be the space of input signals. Defining 

the mapping L : U —> R n as 

Equation 3.1 Lu = // o J <P(t t ,s)B(s)u(s)ds. 

It is obvious that the desired state transfer is possible if and only if the 
equation Lu=d has a solution, i.e. d 6 Im L. 

It is easily verified that L is a linear mapping, but since it does not 
act between two finite-dimensional vector spaces, L does not have a finite- 
dimensional matrix representation. 

Taylor’s ‘Introduction to Functional Analysis’ helps us prove the following 
theorem [Taylor, page 250]. 

Theorem 3.1 If X, Y are complete inner product spaces and L : X — + Y is 
a linear continuous operator then 


Im L = Im LL 
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R n is a Hilbert space but what kind of space is UI Define the inner product 
for U as 

(u,v) u = f u(t) T v(t)dt 
Jtn 


'to 

and it can be proved that U becomes a Hilbert space. The adjoint operator 
L * is determined by 

(Lu, d)nn = d T [ <P(t! , s)B(s)u(s)ds = 

Jt 0 

{5 T (s)# T (^,s)d} u(s)ds = (u, L m d) u 

Vu e u,d e R n 

Consequently we get 

L* : R n -» U as (L*d)(t) = B T {t)<P T (tt,t)d 

and finally 

r‘l 

"o 

We thus have a linear mapping 

LL m : R n R n 


LL’d = 


<P(t I ,s)B(s)B T (s)<P T (t t ,s)ds 

Jtn 


that is given by the symmetric, positive semidefinite n x n matrix. 

W{*o,t,)= r *(ti,s)B(s)B T (s)<P T (t t ,s)ds 

Jto 

Theorem 3.2 We can take a system from x(t 0 ) = x 0 to x(ti) = x t if and 
only if 

d = Xi — $(t , , t 0 )x 0 6 Im W(t 0 , ti ) 

We also have that the control signal u with minimum norm (energy) is given 
by 

u{t) = B T (t)$ T (tj,t)a 
there a is just any solution to 


W(t 0 ,t t )a = d 
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Remark 1 The point with the above is that it is much easier to characterize 
Im W than Im L, because W is an ordinary matrix. 

Proof Let L be as in (3.1) 

(if) Suppose first that d G W(to,t)), i.e. d G Im LL*, then d = Im LL* a 
for some a G R n . Let u = L m a and we get Lu = LL* a = d. Thus, d G Im L 
and the state transfer is possible. 

(only if) Suppose now that the state transfer is possible, i.e. d G Im L. 
Furthermore, suppose that d Im W(to,ti), i.e. d Im LL*. This will 
give a contradiction. 

Recall that for any matrix A it holds that Im A = ( ker A T ) 1 '. Since LL* 
is a symmetric matrix, we get d ^ ( ker LL*) ± . This implies that there is a 
z G ker LL*, i.e. LL*z = 0, such that (z, d) R n ± 0. But, LL*z = 0 implies 
that 0 = (z, LL*z)ru = ( L*z , L*z) u . Hence, L*z = 0. Now the contradiction 
easily follows since 

0 ^ (z, d) R n = (z, Lu) R n = (L*z, u) u = 0 

The final step to prove the optimality of u. Let u be any solution of Lu=d. 
Then Lu = Lu so L(u — u) = 0. This gives 

0 = (a, L(u — u)) R n — ( L* a , u — ti), = («, u — u). 

Hence, (u, u ) = (u, u). We now get by using the Cauchy-Schwartz inequality 
that 

(m, u) = («, u) < («, u) 7 / 2 (u, u ) 1 ^ 2 . 

Dividing by ( u , u) 1 ^ 2 yields that 

(u,u) 7/2 < ( u, u y' s . 


Hence, u is optimal. □ 
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3.1 Reachability for Time- Invariant Systems 

For a time-invariant system 

x = Ax + Bu W{t 0 ,t 1 ) = f 1 e A{tl ~ ,) BB T e AT ^- s) ds 

Jig 

the question about reachability is radically simplified. 

Definition 3.1 Let (A,B) be a matrix pair, where A is nxn.The reachability 
matrix V is defined as 

[B,AB,...,A n ~ 1 B] 

Theorem 3.3 For all (to,^i) such that to < ti we have 

Im W{t 0 ,ti ) — Im r 

Proof I. Im r C Im W 

Im TCIrn W & ( ker r T ) x C ( ker W T ) ± & ker W x C ker 

Presume that a € ker W , i.e. Wa=0 so a T Wa = 0 and hence it follows that 

a T e M t i-A B = 0 Vs 6 [to,lt]> 

Derivation with regard to s a couple of times and s := t t gives 
a T B = 0 . . . a T AB = 0 a T A n ~‘ B = 0 i.e. a (E ker r T . 

II. Im W C Im T 

In the same way as above we are going to prove that ker T T C ker W. 
Suppose that a T F = 0. By Cayley-Hamilton follows 

a T A k B = 0 k = 0,1,2,... 

Accordingly we have 



So it follows that a T W — 0, i.e. Wa=0 , i.e. a 6 ker W. □ 

Remark 2 Since Im T = Im VF(<o,^i) for any interval (f 0 , *i), we see that 
in the time-invariant case the image of the reachability Gramian is inde- 
pendent of the interval (to,ti). However, this does not imply that the state 
transfer can occur during a fortuitous short time interval. 
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Definition 3.2 Let n be the dimension of the state space. The pair (A,B) 
is said to be completely reachable ifT has full rank, i.e. 

rank T = n 

Definition 3.3 The reachable subspace 72 is defined as 
H = Im [B,AB,A s B,...,A n ~ 1 b] 

We easily see that 72 is the set of states that can be reached from the origin. 
Lemma 3.1 The reachable subspace 72 is A-invariant, i.e. 

Alien 

In particular, e M n C n for all t £ RA . 

Proof Since, by the Cayley-Hamilton theorem, A n is a linear combination 
of A * for j—0,1,. . . ,n-l it follows that 

AK= [AB,A s B,...,A n B\ C Im [fl, AB , . . . , A n ~‘ B] = K. 

Moreover, by induction we get A*H C n, which implies that 

e At n = £ ^A j n c n □ 

j=o 

To further clarify the picture we note that if the state of the system is in n , 
at some instant, it is impossible to steer the state out of n. Neither is it 
possible to enter 72. from an initial state not in 72. Particularly we have that 
if x 0 , x t € 72 then the state transfer can occur in just any time t. The points 
that can be reached in a time t from a given x 0 establish the plane 

n(x 0 , t) = e M x 0 + 72. 
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4 Stability 

A very essential problem when designing a control system is how to avoid 
instability, i.e. that the output increases without limit. The following sec- 
tion is an abridgement of chapter four in “An Introduction to Mathematical 
Systems Theory” by A. Lindquist and J. Sand [Lindquist/Sand]. All theory 
dealing with the alternative approach, the Lyapunov equation, is omitted. 

Intuitively an input-output system is stable if a bounded input produces 
a bounded output or if the output tends to zero, or at least remains bounded, 
when the input is zero. 

For nonlinear systems, stability in this sense is typically dependent on the 
initial conditions and the specific input applied. Hence, in general, stability 
is not the property of a system, but rather the property of a solution. 

This chapter deals only with the stability of time-invariant linear systems, 
a subject which is drastically simplified by the fact that the complete set of 
solutions of the system x = Ax can be displayed explicitly by means of the 
Jordan form. As a consequence, it is enough to check the eigenvalues of A 
in order to determine whether a bounded input produces a bounded output, 
and thus it will be meaningful to talk about stable systems. 

4.1 Stability of Continuous-Time Systems 

We want a bounded input to give a bounded output, which is sometimes 
abbreviated as H/HO-stability. 

Definition 4.1 The system 

( x(t) = A(t)x(t ) + B(t)u(t) 

{ y(t) = C(t)x(t) 

is input-output stable if there is a k such that 

||u(°<)||<l <G [<o,oo) } ^ * € [*0’°°) 


for every t 0 . 
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Example 4.1 Consider the time-invariant case, where (A,B,C) are constant 
matrices. Then 

y(t) = r Ce A ^Bu{s)ds. 

J to 

Defining G(t) = Ce At B, 

l|y(OII < / ||G(f-.s)|| ||u(s)||ds ( since ||u(f)|| < 1) 

Jto 

< f" \\G(v)\\d<, 

J *0 

< l|C||||B|| f~'° He^lK, 

Jto 

i.e., a sufficient condition for input-output stability is that the integral 
Io° ll e ' 4< || c ^ 15 convergent. □ 

4.2 Stability matrices 


Let us study the homogeneous system: 

Equation 4.1 x = Ax; x(0) = x 0 . 

Definition 4.2 The system (f.l) is stable if the solution is bounded on the 
interval [0,oo) for all initial values x 0 and asymptotically stable if x(t) — > 0 
when t — ► oo for all xq. 

Theorem 4.1 (1) The system (4-1) is asymptotically stable if and only if 
the real parts of all the eigenvalues of A are less than zero, i.e. the eigenvalues 
are all located in the open left half plane. 

(2) The system (4-1) is unstable if A has at least one eigenvalue in the open 
right half plane. 

Proof In this proof we shall use a fundamental result from linear algebra, 
the Jordan decomposition theorem. This theorem guarantees the existence 
of a basis for 72." in which the representation of the linear mapping A takes 
a particularly simple form. 
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Transform the matrix A to Jordan form A = TJT 1 , where J is a block- 
diagonal matrix. 

J — diag{J i , Jj, . . . , J r ) 
and each d v x d„ block J v has the form 


Jv = 


A„ 1 
A v 


0 


A„ being an eigenvalue of A. Thus, 


r e Jit 


e At = T 


oJit 


0 ‘ 

1 ’ 

A„. 

0 1 


so it remains to analyze each e Jvt . But J v has the form 

J v = A „/ + S v 

where S v is a shift matrix 


S v = 


0 1 0 

0 1 

o 

•. 1 

0 0 


of dimension d v x d v , having the property that S' = 0 for i > d v . Conse- 
quently, 


e Jvt _ e *vt e s v t _ e A •< | / + f5 + l 7 5 4 + ... + 7 -i— —S iv ~ 1 } 


t* 

Jr 


(d v 


and therefore, setting o v — Re A* and ui v = Im \ v , 



4 STABILITY 


13 


Equation 4.2 e Ai = £„ e" vt P v (t)(cosu v t + sinu v t), 

where P„(i) is a matrix- valued polynomial of dimension d v — 1 in t. From 
this expression it follows that (1) e At x o * 0 for all x o if and only if 
a v = Re\ v < 0 for all v and that (2) e At x 0 -» oo for at least one x 0 if some 

a v > 0. □ 

Lemma 4.1 The system in equation 4.1 is stable if and only if all eigenvalues 
of A are located in the closed left half plane and any eigenvalues on the 
imaginary axis correspond to one dimensional Jordan blocks. 

Proof By theorem 4.1 (1) we only need to worry about terms in (4-2) for 
which cr v = 0, i.e. e° vt = 1. These terms will remain bounded if and only if 
the degree of P v is zero, i.e. d v — 1. 

Definition 4.3 A is a stability matrix if Re A(A) < 0. 

Theorem 4.2 If A is a stability matrix then the time invariant system 

{ x = Ax + Bu 

y = Cx 


is input-output stable. 

Proof If all eigenvalues of A have negative real parts so that all cr v in (4-2) 
are negative then 

r°° 

/ \\e At \\dt <oo 

Jo 

and hence, in view of example 4.1 the system is input-output stable. □ 
The last theorem is very important for us because it deals with the kind of 
system we use when modeling an aircraft. 
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5 The Systems 

The test was accomplished with three systems of different kinds. All systems 
use a single input and produce a single output, a so called SISO- system. 
As we will see later all systems have the necessary property of complete 
reachability, as was discussed in chapter 3. 

That a bounded input gives a bounded output is sometimes abbreviated 
as BIBO- stability. This highly desirable property for a system was discussed 
in chapter 4 and will be further examined for each specific case. 

5.1 Transfer Functions 


This subject is discussed in Etkin’s book “Dynamics of Atmospheric Flight” 
[Etkin, page 50-51]. He writes 

System analysis frequently reduces to the calculation of system 
outputs for given inputs. A convenient and powerful tool in such 
analysis is the transfer function, a function G(s) of the Laplace 
transform variable s [Complex valued], that relates input u(t) and 
output y(t) as follows, 


G{s) = y(s)/u(s) 

where ( ) denotes the Laplace transform. So long as u(t) and 
y(t) are Laplace transformable the transfer function defined above 
exists. However, it will in general be a function of the initial 
values of y and its derivatives, and moreover, for nonlinear and 
time varying systems, of the particular input u(t) as well. Such 
a transfer function is of relatively little use. We can however 
obtain a unique function G( s) if (I) the system is linear and time 
invariant, and (II) it is initially quiescent, i.e. at rest at the origin 
in state space with no inputs. 

He continues, 


When u(t) and y(t) are zero for t < 0, the Laplace and Fourier 
transforms are simply related, i.e. u(iu) = U(u). It follows that 

Y{») 


G(iu) = 


u M 
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Sometimes it is G(iu > ) that is called the transfer function. 

When examining different transfer functions in a Bode diagram it shows 
that there are systems with the same absolute value curve but with different 
phase curves. Of all systems with the same absolute value curve there is 
one with less negative phase advance, it is called a minimum phase system 
[Glad/Ljung, page 109]. 

I give the following theorem without a proof. 

Theorem 5.1/1 theorem with a rational transfer function is in minimum 
phase if and only if it has neither poles nor zeros in the open right half plane. 

The others are called non minimum phase systems. This distinction is very 
important because we know from one dimensional control theory that a sys- 
tem with zeros in the numerator will start off in the opposite direction. This 
bad quality can make the system difficult to control. 

A A-value less or equal to zero are assumed in the following calculations. 


5.2 System 1 


i = 


0 1 

0 A 


x + 




gives the system dynamics y = Ay + u 
The reachability matrix 



has full rank for all lambda and the system is therefore completely reachable. 
System l’s transfer function 


Y(s) 


1 

s(s — A) 


U(s) 


without neither poles nor zeros in the open right half plane indicates that it 
is a minimum phase system and should therefore be easy to control. 
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The Jordan transform of the matrix A is P 1 JP where 


o 

0 

1 

p — 

A/(A — 1) 1/(1 -A)' 

o 

V- 
1 

7 r — 

0 1 


Because all eigenvalues of A are located in the closed left half plane and the 
eigenvalue on the imaginary axis correspond to an one dimensional Jordan 
block we know that the system is stable. Owing to this quality we are guar- 
anteed that when a fortuitous in signal ultimately equals zero, the solution 
to the system x = Ax is bounded on the interval [0,oo). This of course im- 
plicates that also the output is bounded. Referring to previous theory the 
eigenvalue on the imaginary axis prevents input-output stability. 


5.3 System 2 


‘ 0 

1 

0 

0 ■ 


' 0 ' 





y 

\ 

0 

A1 

1 

0 

X + 

e 

w 

II 

o 

o 

o 

X — 


y 


0 

0 

0 

1 

0 



u 


. 0 

0 

0 

A2 . 


. 1 . 





u 

/ 


gives the system dynamics y = \1 y + u -f ew u — \2u + w. 
The reachability matrix 


’ 0 
e 
0 
1 


e eAl eAl 2 + 1 

eAl eAl 2 + 1 eA 3 + A1 + A2 
1 A2 A2 2 

A2 A2 2 A2 3 


has full rank for all values on A and e and the system is therefore completely 


reachable. 

System 2’s transfer function 


Y(s) 


es 2 — eA 2s + 1 . 

s 2 (s — A l)(s — A 2) 


gives for negative A’s that there are no poles in the open right half plane. 
The numerator es s - eA 2s + 1 = 0 give the solution 


A 2 , /A 2* T 

, -T ± VT“T 
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L h Af^r v ^K n r* atiVe 1 th ? "" haVe 3 2610 iD ‘ he °P en ri S ht half Plane. 

left half nlane V 7 "“"“l' ‘° the P °' eS f ° r a P™*™ 3 are in the 
eft half plane. So the system should be easy to control for a positive e and 

p o a y more difficult for a negative value on the variable. 

taking the Jordan transform of A gives that J equals 

T A2 0 00 

0 Al 0 0 

0 0 0 1 

0 0 0 0 

Crying through the same discussion as for system 1 we see that system 2 
a* the same properties, stable but not guaranteed input-output stable. 

5.4 System 3 


X = 


' 0 

1 

0 

0 

0 


' o ‘ 




0 

Al 

1 

0 

0 


e 



f y ^ 

0 

0 

0 

0 

0 

0 

1 

0 

0 

1 

X + 

0 

' 0 

w y = 

1 o 0 0 0 ] X X = 

y 

U 

0 

0 

0 

0 

A2 


1 



U 










\ / 


Kxy UCU1II 

The reachability matrix 


«= X2u + w 


f = 


0 

e 

0 

0 

1 


6 

eAl 

eAl 2 

eAl 3 + 1 

eAl 

eAl 2 

eAl 3 -f 1 

eA 4 + Al-f A2 

0 

1 

A2 

A2 2 

1 

A2 

A2 2 

A2 3 

A2 

A2 2 

A2 3 

A2 4 


has always full rank and the system is therefore completely reachable. 
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System 3 s transfer function 


Y(s) 


- t\ 2s* + 1 

's J {? - (17 + A 2)s + Tl • A J) 


can be examined by Routh’s algorithm, 
gives the following tableau. 


The numerator es 3 — e\2s s + 1 


e 0 ' 
— eA2 1 
1/A2 0 

L 1 0 


e e > ;• °: € t> ? ;; :r\ > „ a r ,ive ■ that the wt 

“r the rr has r — ^he 0 OP t“s.irJXT r 

a negltive £ *- that the ^ h “ 

The solution to the denominator 


- 


+ AZ)S + XI - A 2) = 


gives that for all negative values on \1 and \2 we have three zeros on the 
aginary axis and two zeros in the open left half plane. If either \1 or \2 
equals zero we get four zeros on the imaginary axis and one in the open left 

All this together gives that system 3 never will be a minimum phase 
system and will therefore be more difficult to control. ? 

iaking the Jordan transform of A gives that J equals 


A2 0 0 0 0 

0 A1 0 0 0 

0 0 0 1 0 

0 0 0 0 1 

0 0 0 0 0 


and this implies that system 3 is stable. The always present eigenvalues that 
equal zero prevent the system to be guaranteed input-output ftable. 
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6 Derivations 

The fundamental idea of the test was to simulate an aircraft that is directed 
by the flight control to follow a certain path. The given trajectory can be 
seen as a set of points that shall be passed at a certain time. By way of 
introduction our intention was to determine how exact a given path in three 
dimensions could be tracked with maintained comfort for the passengers. 
Priority one was to minimize the acceleration and thereby the stress to the 
individuals. 

To start with, the test was accomplished in one dimension. This sim- 
plify the calculations radically and is a common approach to such experi- 
ments. Given the set of points {yo, yi , • • • > 2/n} an d the corresponding time 
{ t 0 , t t , . . . , t n } , we would like to perceive the control laws { u 0 , u t t } 

that take the system through the points in such a pleasant way as possible. 

Consider the control u k that takes the system from state vector x k to 



Because t E [t k , t k +i ] the state of the system will be 

z(*) = e Mt ~ tk) x k + [ e A{t ~ s) Bu k (s)ds 
Jit 

and as the state of the system is x k+l at time * i+ , we receive the condition, 

Equation 6.1 

Xk+i = e A(ik +‘~ tk) x k + f t+1 e A ( tk+1 ~^Bu k (s)ds. 


The solution u k to equation (6.1) that minimizes the energy of the control 
signal is then given by, see chapter 3. 

Equation 6.2 

u k (t) = B T e' ATi ^ <t+i e~ Aa BB T e~ AT, ds ^ x k+1 - e~ Mk x k ) 
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stapiZll a b ; v ;:; convenient if the **** - <«*> - uW „« 

of See OUr path “ the aircraft is flo"n through a largo number 

of pomta by an autopilot. With a specific time interval the plane Reives a 

The time intevai /“ ^ VehkIe ^ ““ Path with hi « h a <™racy. 
the automatic pilot works w“h' C °" Stant ^ determine<1 * tha fluency 

Assumption 6.1 Let t k+l -t k = h. 

Assumption (6.1) can be used to simplify the integral in equation (6.2). 

- 4 


L 


*k+l 


S BB' 


*ds = {t = s- t k ] = 


"-At, 


J 0 

v 


Definition 6.1 


matrix constant 


M = [ e~ Ar BB T e~ ATr dr 

Jo 

The equation (6.2) can be rewritten as 

Equation 6.3 


“*W = i 4+ , - I( ) 


a kn s rtio w :. have to app,y — kipd o f conditCrr^^^r 

The control it is the actual control that the pilot or the auto-pilot achieve 
ery natural choice is therefore to require that the control a is continuous' 
Assumption 6.2 


Mh+i) = u k+I (t k+I ) 
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By this we acquire (n-1) conditions and allow us to write 

«*(£*+/) = «i + /(<t +J ) 

This equation can be simplified by 

Definition 6.2 

Z = M~* e~ Ah 

W 4 e- ATh M- l e~ Ah + M~ l 

and finally we obtain the modified expression 

B ( Z T xt — Wxk+i + Zxit+z) = 0 k = 0, 1 , . . . , n — 2. 

Written in block diagonal form it becomes, 

Equation 6.4 

n ~ XO ■ 

Z T -W Z ... 0 0 0 1 X, 

0 -W ... o 0 0 x 2 

B T ; ! I : : : : =0 

o 0 0 ... —W Z 0 x n _ 2 

.° o 0 ... Z T -w z\ x n _! 

. x n . 

As our three systems are very different the calculations will differ from here 
and all further computations have to be treated separately. 
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6.1 System 1 


This is the original system that was first implemented in Matlab. Each state 
vector x k consists of two parts, a known coordinate y k and an unknown ve- 
locity y k . By partitioning the matrixes in definition 6.2 as 

Z = 211 212 Z T 

z 2l z 22 ’ 

and using the notations given in the following definition, the unknowns can 
be kept on the left hand side and the given position coordinates can be moved 
to the right hand side. 

Definition 6.3 


*11 

*21 

, w = 

w u 

w 12 

*12 

*22 _ 


w 21 

W 22 


W, B = B T 

w 21 
w 22 

= [wez] 

= B T 

z 12 
222 

= M 

z$ = B T 

z 2\ 

z 22 

= M 

W r B = B t 

\ ^11 
Wn 

= [wti] 

Z b = B t 

Z 11 
*21 

= M 

z* = b t 

*11 

*12 

= M 
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We get the system 







_ [ 


yo 


z?i - 

-W, B 

Zl . 

0 

0 

0 



y\ 


0 

zft 

-W B 

0 

0 

0 



y 2 


0 

0 

0 

.. -w? 

k 

0 


J/n- 2 


0 

0 

0 

7 b 

A ll 

-W B 

Zt. 


J/n— 1 









y n 









yo 

- 

■ -Zrl 

w t b 

-Z * » 

0 

0 

0 



y i 


0 

-zf x 

W B 

0 

0 

0 


y 2 


0 

0 

0 

... W B 

-Z* 

0 


yn-'. 

2 

0 

0 

0 

... -Z$ 

W B 

-Z? u . 


Un - 1 








L y n 

J 


The right hand side consists of known parameters and is therefore a constant 
vector. Our system has (n-t-1) unknowns but only (n-1) equations so we need 
two more constraints. 

After reading Per Enquist’s paper “Control Theory and Splines, applied 
to Signature Storage” [Enquist] I decided to use the natural boundary con- 
ditions, y 0 = 0 and y n = 0. This can be seen as a very real behavior for a 
vehicle and has also given the best results in former experiments. Enquist 
writes “This will let the initial direction and constant velocity of the system 
be decided so that the control energy is minimized” [Enquist, page 16-17]. 
The system dynamics equation y = Ay + u gives 

A y 0 + u 0 = 0 where u 0 = B T M~‘ (e^x, - x 0 ) = B T Zx t - B T M~ l x 0 . 


U B = B T 


-i 

-1 

mn 

m 12 

-i 

-1 

m 21 




Definition 6.4 
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We get 

(A - u*) y 0 + z*y t = U B y 0 - Z B y t 

The system dynamics equation together with earlier definitions give 

+ u n ..,(t n ) — 0 where u n _/ (f„) — B T e~* T h M~ l [e~ Ak x n — ar„_;) = 

B T e~ ATk Zx n - B T e- ATk M~ 1 x n . 1 = B T Wx n - B r M~ l x n - B T Z T x n . t 

We get 

~Zu y n -i + (A + W[ B — Ug) y n = Z B y n - t + (uj^ — y n 

Add these two equations to our considered system and the number of un- 
known parameters equals the amount of equations. Thus is the problem 
solvable and the Matlab program that uses Gaussian elimination is displayed 
in mprl21der0.m. 

6.2 System 2 


By this system a new approach was introduced for the convenience of the 
passengers. Instead of direct using the performed control signalu we use its 
derivatives to control the aircraft. The formula ii = X2 u -|- w gives the con- 
nection between the control u induced by the pilot and the artificial control 
w that actually flies the plane. 

Each state vector xt consists of the known coordinate yt and the unknown 
parameters, velocity yi, control signal 14 and its first derivative tit. As we 
have (n+1) state vectors and each state vector consists of three unknowns 
it becomes a total of 3(n-|-l) unknown variables. The constraint that we 
require the control signal u to be continuous gives only (n-1) conditions. If 
the restrictions are introduced that also ti and u have to be continuous, we 
get further 2(n-l) conditions. 

Assumption 6.3 


Uk{tk+l ) = Uk+i ( tk + i ) 
Uk{tk+l) = U k+I (t t+1 ) 
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Applying this to equation (6.3) becomes for the first condition 

B t A t ( Z T x k — Wxk+t + Zxk+ 2 ) = 0 k = Oy 1 . yti — £ 
and for the second condition 

B t A t A t (z T x k - fVxjt+j + Zz i+ ,) = 0 k = 0, l , . • . , n - 2. 

By partitioning the matrices in definition 6.2 as 

Z\\ Z 21 z 3l z 4l 

z \2 z 22 z 32 z 42 

Z 13 z 23 z 33 z 43 

. z \4 z 24 z 34 z 44 . 

W 14 " 

W24 
W34 
W44 . 

and using the notations given in the following definition, the unknowns can be 
kept on the left hand side and the given position coordinates can be moved 
to the right hand side. The matrix notation . symbolizes a whole row or 
column. 

Definition 6.5 

Continuous control signal, u: 

W? = B t [w.s w.s w. 4 ] = [ew 22 + w 42 ew 2S + w 43 ew^ + w 44 ] 

Z,* = B t [z, 2 Z'S z. 4 \ = [tz 22 + z 4 i ez 23 + z 43 cz 24 + z 44 ] 

Z,f = B t [z s . z s . z 4 ] = \tz 2i + z 24 tz 32 + z S4 tz 4S + z 44 ] 

W r B = B t [wj] = \ew 21 + w 41 ] 

Z = B t [z a \ = [ez 21 + z 4] ] 

Zfj = B T [z,] = [ez t2 + z 14 ] 


Z = 


z ll 

z 12 

z 13 

z 14 

Z 21 

222 

z 23 

z 24 

Z 31 

^32 

z 33 

z 34 

. z 4l 

2 42 

z 43 

z 44 . 


Z 1 = 


w = 


W U W 12 w 13 

W 2 1 W22 w 23 

W 3 1 W 32 w 33 

L W41 W42 w 43 
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Continuous first derivative of control signal, u: 

W* = B t A t [w. £ w.s w. 4 ] = [ew JS + eXlwtt + w 3S + X2w 4t 
ew ‘s + tXlwgs + w 33 -f X2w 43 ew, 4 + eXlw t4 4- w S4 + A 2w 44 ) 

= B A [z .2 z ,3 z, 4 ] — [ez/g + eXlzgg + z 3 g -|- X2z 4 g 
€Z is + tXlz i3 4- z 33 4 A 2z 4 s tz 14 + tXlzg 4 + z 34 + X 2 z 44 ] 

Zu=B t A t [z s . z 3 , z 4 .] = [ez gl + eXlzgg + zg 3 + X2 z S 4 
ez si + eXlz ss + z 3S + A 2z S4 ez 4 i + eXlz 4S + z 43 + A 2z 44 \ 


= B A [w»./] = [ewji + eXlwgt 4- Wst + X2w 44 ] 
^r« BA [zj] = [ezij 4- tXlzgi -|- z 3 i 4* A 2z 4 j] 

^ rt B A [zj ,] = [ezjj 4 eXlzjg 4- zj 3 4* A 2zj 4 J 

Continuous second derivative of control signal, u: 

W t B = B T A T A T [w.2 W.3 W. 4 ] = 

[eXlw I2 + (eA I s 4 l)w 22 4 A 2w 3S + X2 s w 4 g 
eXlw i3 + (tXl 2 + l)w S3 + X2w 33 + X2 s w 43 
eXlwj 4 + (eA I s + l)wg 4 -f A 2w 34 +A 2 s w 44 ] 

Z* = B T A T A T [zg Z' 3 z 4 ] = 

[tXlztg + (eA I s + 1 )z S g + A 2z ss 4- X2 s z 4 g 
eXlzis + (eAi + l)z 23 4 A 2z 33 + A 2 s z 43 
eXlz i4 + (eA I s 4- 1 )z t4 4- A 2z 34 +X2*z 44 ] 
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Zfi = B t A t A t [z s . z 3 . z 4 .] = 
[eXlzsi + ( eXl 2 + 1 )z 22 + AUz 2 j + X2 S z 24 
t-Xlzgi “t* {(-Xl -j- l)z 32 “f - X2z 33 + A 2 2 z 34 
tXlz 41 -f (eA I s -f 1 )z 42 4- X2z 43 + Ai? 2 z^] 


—B A T A t [w,i] — [eXlwu + (cXl e + l)w 2J + X2w 3t + Ai? 2 w 4 t\ 

Z* = B t A t A t [z. i ] = [e\lz n + (eA7 2 + l)z 21 + X2z si + A 2 s z 4t ) 

Zf\ = B t A t A t [z ! ] = [eAiz„ +{eXl 2 + 1)z is + X2z 13 + X2 s z I4 ] 

We get the following systems, written in block diagonal form. 

Each state vector is divided into two parts, a known portion x£ which 
contains the given position y* and an unknown portion x% that contains the 
parameters y±, u * and u*. All right sides consist of known variables and are 
therefore constant vectors. 

Continuous control signal it: 


' z B 

-W B 

Zl ... 

0 

0 

0 



X V Q 

x v . 

0 

Z$ 

-w B ... 

0 

0 

0 


x v 2 

0 

0 

0 

-W t B 

Zl 

0 


X n- 2 

0 

0 

0 

Zu 

-W B 

Zl 


X n-1 








K J 







X p 0 

'-Z?, 

W r B 

-Z t b u ... 

0 

0 

0 

1 


x\ 

0 

-z?i 

w r B ... 

0 

0 

0 


x p 

0 

0 

0 

W B 

-Z* 

0 


T P 

x n-2 

0 

0 

0 

-Z*, 

W B 

-Z* 

J 


T P 

1 
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Continuous first derivative of control signal, ti: 



Continuous second derivative of control signal, u: 




6 DERIVATIONS 


29 


Having totally 3(n+l) unknowns but only 3(n-l) constraints we chose to 
enter differential approximations for the first and last state vector, this will 
decrease the number of unknown parameters by six and thus make the prob- 
lem solvable. Remember that the time interval — ft is constant and 
represented below as h. 

The needed velocities are approximated as: 


yo 


yi - yo 


h 


Vn ~ yn-t 

Vn = h 

The needed control signals axe approximated as: 



Vn-t = 


Uo = 


u n = 


i/n-/ ~ yn -2 

h 

yo - y/ 

h 

ifn—1 Vn 


The needed first derivative of the control signals are approximated as: 

Vs ~ yi 

yz - — 7 — 


y n -z = 


U] = 


y n -2 - Vn-s 

h 

yi - yz 


Un-l = 
Uo = 
U n = 


yn-Z - yn -/ 

h 

U 0 ~ U t 


U n -1 ~ U n 


The Matlab program that solves the task for system 2 using Gaussian elim- 
ination is displayed in mprl41knovel.m. 



6 DERIVATIONS 


30 


6.3 System 3 


This system uses the same approach for the convenience of the passengers 
as system 2 does. The only difference is that we now also use the third 
derivative of the control signal u to actually fly the plane. The formula 
u== \2ii ' + w gives the connection between the control u induced by the pilot 
and the artificial control w that actually flies the plane. 

Each state vector consists of the known coordinate yu and the unknown 
parameters, velocity y*, control signal u*, its first derivative and its second 
derivative it*. As we have (n-f 1) state vectors and each state vector consists 
of four unknowns it becomes a total of 4(n+l) unknown variables. The con- 
straint that we require the control signal u to be continuous gives only (n-1) 
conditions. If the restrictions are introduced that also ri, u and u have to be 
continuous, we get further 3(n-l) conditions. 

Assumption 6.4 

Uk(tk+i ) = Uk+l ) 

Uk(tk + l ) “ tyfe+i (^Jfc+i ) 

(tk+i ) =£*+i (**+/) 

Applying this to equation (6.3) becomes for the first condition 

b t a t ( z T Xk _ Wxk+1 + Zxk+t) = 0 k = 0, 1 , . . . ,n — 2, 

for the second condition 

B T A T A T (Z T x t - Wx k+i +Zx t+l ) = 0 k = 0,l,...,n- 2 
and for the third condition 

B t A t A t A t (z T x t - Wx k+I + Zx i+S ) = 0 k = 0, 1 . .,n — 2. 

By partitioning the matrices in definition 6.2 as 


Zn 

Z\2 

Z 13 

Z 14 

Z\5 


z 11 

Z2\ 

Z31 

Z41 

Z51 

Z 21 

z 22 

z 23 

Z 24 

z 25 


■212 

Z22 

Z32 

Z42 

Z52 

Z 31 

Z32 

Z 33 

Z 34 

Z 35 

Z T = 

■213 

Z23 

Z 33 

Z43 

Z 53 

Z 41 

Z 42 

Z 43 

Z 44 

Z 45 


-214 

Z24 

Z34 

Z 44 

Z54 

z Sl 

Z 52 

Z53 

Z54 

Z55 


. Z l 5 

Z25 

Z 35 

Z 45 

Z 55 
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W 11 

W U 

™13 

™X4 

^15 

W 21 

W 22 

™23 

™24 

W 2 5 

^31 

W32 

^’33 

™34 

™35 

™41 

W42 

™43 

™44 

™45 

™51 

^52 

^53 

™54 

W5S 


and using the notations given in the following definition, the unknowns can be 
kept on the left hand side and the given position coordinates can be moved to 
the right hand side. The matrix notation . symbolizes a whole row or column. 


Definition 6.6 

Continuous control signal, u: 

W[ B = B t [w,2 W. 3 w 4 w. 5 ] = 

[ew 22 + ru'ss tw zs + w ss ew £4 + w s 4 tw ss + w jj] 

Z B = B t [z . 2 z.s z .4 z. 5 ] = [tz 22 + z 5S ez 23 + z 53 tz S4 + z 54 ez 25 + z 55 ] 
Z B = B t [z 2 . z 3 . z 4 . z 5 ]-[ez 22 + z 25 ez 32 + z 35 ez 42 + z 45 ez 52 + z 55 ] 

W r B = B t [wj] = [ ew st + wsi] 

Z B U = B t [z.i] = [ez 2t + zsi ] 

Z B = B t [zj ] = [ez 12 + Zjs] 

Continuous first derivative of control signal, u: 

W B = B t A t [w. 2 w.s w, 4 w. 5 ] = [ ew t2 + e\lw 22 + w 42 + X2w 52 
ew 13 + e\lw 23 + w 43 + \2w 33 ew 14 + e\lw 24 + w 44 + X2w 34 ] 


Z B = B t A t [z, 2 z .3 z. 4 z.s ] = [tz 12 + eXlz 22 + z 42 + X2z 52 
ez 13 + eXlz 23 4- z 43 + X2z 53 ez i4 + eXlz 24 + z 44 + X2z 54 ] 

Z B = B t A t [z 2 . z 3 . z 4 . zs] = [ ez 24 + eXlz 22 + z 24 + A 2z 23 
tz 3l + tXlz 32 + z 34 + X2z 35 ez 4 i + eXlz 42 + z 44 + X2z 4S \ 
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W r — B — [cwu A cXlwsi 4" + A 2wsi] 

= B t A t [zj] = [ezjj 4- e\lz2i + z 4 i + A 2zsi] 

Z T i = B^ A T \zt\ = \tzu 4- e\lzt2 4- Zjj + A 2zi$\ 

Continuous second derivative of control signal , ii: 

W B = B t A t A t [W' S W'S W.4 w.5] = 

\e\lW 12 4“ 6 A I s W 22 + 4“ \ 2 wj 2 “h A 2* W $2 

eXlwjs 4~ e\l s ti)23 4- 4* A 2w^s 4~ A 2*w$s 

6 A lwi£ 4" fA 1 W24 4- W34 4“ A 2w^4 4" Ai?^it?5^] 

Z£ = * r i4 r j4 r [z. f z, z m4 *.,] = 

\t\lz12 4" 6 A l 2 Z22 4- ^5 4~ X2z^2 4“ A 2* z$2 
6 A 1 Z \3 4- eAi zgj 4" 2*5 + X 2 Z 43 4* X2* Z 53 
e\lzi 4 4" cAi z^ 4“ 4~ X2z 44 4~ A 2 s z$ 4 ] 

Zfi = 5 r A r v4 r [z 5 . z*. z^. z 5 J = 

[eAiz^i 4- tXl 2 Z22 4- 4- Ai?z^ 4- AiJ^z^j 

eXlzsi 4 " tXl^ Z32 + Z33 4 “ X2Z34 4* A 2* Z35 

tXlz^i 4" cXl Z42 4~ Z 43 “1" Ai?z^ 4" Aj&^z^] 

= 2J r i4 r j4 T [u7j] = [tXlw u + tX I s w 2 i 4- W31 4- X2 w 41 4- X2 s w 5J ] 
^r% = B T A T A T [z,i ] = [t\lz u 4“ eXl 2 Z21 4* zsj 4* A 2 z 4 i 4- X2*z$i] 

Z r i — & A A [zj J = [eXlz u 4" tXl 2 Z 12 4" zis 4- A 2zj 4 4“ A 2 s Zj^J 
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Continuous third derivative of control signal , u: 

B t A t A t A t [w.s w.s w.i w.s] = 

[eX l 2 w ts + ( eXl 3 A 1 )w ss A X 2 w 3Z 4- A 2 2 w 42 A X2 S w 5S 
eXl s wis + (eXl 3 A 1 )w 2 s + X2w$s A X2 2 w 4 s A A 2 W 53 
tXl 2 Wu + (eXl 3 A 1 )w 2 4 4- A 2wsj + A 2 w 44 + A 2 wn\ 

Zi= B t A t A t A t [z .2 z.s U 2.5] = 

[eA i 2 Z 12 A (cA 1 3 + 1 ^)Z 22 4- A 2 z $2 A X2 z 4Z A X 2 z$ z 
eX I s Zis + (eXl 3 + 1 )z Z s + X 2 z $3 + A 2 z 4 $ + A 2 z$s 
eX l 2 zj 4 -}■ (cA 1 3 + 1 )z 24 4" X 2 zs 4 4" A 2 z 44 A A 2 z s 4 ] 

zf,= B t A t A t A t [z2. z s . z 4 . z s ] = 

[eX 1 s Z 21 4" (cA l 3 4" 1 )z 22 4" X 2 z 2 s 4" A 2 z Z4 4" A 2 Z 25 
eX 1 2 Z 31 4- (cA I s 4* l)z 32 4” A 2zss 4~ A 2 zs 4 4" A 2 zs$ 
eXl 2 z 4I 4" (cA 1 3 4" l)z 4 2 4* X2z 4 s A X 2 z 44 A X2 z 4 s] 

M>f= B t A t A t A T [w.t} = 

[cA I s Wu 4- (eXl 3 A 1 )w 2 t 4- X2wsi 4- X2 2 w 4 t A A 2 wsi\ 


Z B n =B T A T A T A T [z,} = 

[eXl 2 zu + (cA I s 4- 1 )z 2 i A X2zsi 4- X2 S z 4 i 4* A 2 zsi\ 
Zri= B T A T A T A T [z t ] = 

\tXl 2 Zn A (eXl 3 A 1 )zj 2 A X2zjs 4" A 2 z 44 A X 2 z/j] 
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We get the following systems, written in block diagonal form. 

Each state vector is divided into two parts, a known portion x* which 
contains the given position t/i and an unknown portion that contains the 
parameters t/t, u*, ti* and ti*. All right sides consist of known variables and 
axe therefore constant vectors. 

Continuous control signal u : 








r 

< 1 


yB 

“ll 

-W, B 

Zl 

0 

0 

0 


x\ 


0 

Zft 

-W B 

0 

0 

0 


X 2 


0 

0 

0 

... -W, B 

yB 

0 


<-2 


0 

0 

0 

... zg 

-W B 

z? u _ 


X n- 1 








L 

*n j 








r p 

-1 

\~Z B , 

w B 

-Z* 

0 

0 

0 

1 


T P 

x \ 


0 

-Zf { 

wf 

0 

0 

0 


x\ 


0 

0 

0 

... W r B 

-Z* 

0 


x n- 2 


0 

0 

0 

... -Z B 

W r B 

-Z B 



T P 

X n- 1 








x p 
L •*'n 

m 
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Continuous first derivative of control signal, u: 


\za 

-W B 

Zl B u •• 

0 

0 

0 

0 

zp 

-W B .. 

0 

0 

0 

0 

0 

0 

. -W B 

zP 

0 

0 

0 

0 

• Z% 

-w f 

7 b 

^lu J 


' ~Z r B l 

wf 

-Zf u 

0 

0 

0 

0 

-Z? t 

wf 

0 

0 

0 

0 

0 

0 

... Wf 
... -Z B 
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Continuous third derivative of control signal, u: 
r ...b B 

Zu -W, 
o Zu 

0 0 

0 0 

r ...B ■■■ B 

Z rt W r 

o -Z* 

0 0 

0 0 

Having totally 4(n+l) unknowns but only 4(n-l) constraints we chose to 
enter differential approximations for the first and last state vector, this will 
decrease the number of unknown parameters by eight and thus make the 
problem solvable. Remember that the time interval tk+i — it is constant and 
represented below as h. 

The needed velocities are approximated as: 

yt ~ Vo 
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Vn - J/n-/ 

y " = h 

The needed control signals are approximated as: 


Vi 



y n — i 


y n „i - y n - 2 
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u 0 = 
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u n = 


Vn — 1 J/n 


The needed first derivative of the control signals is approximated as 

yt ~ Vi 


ys = 


Vn-t = 


u, = 


Vn — t Vn — 3 

h 

i/i - v* 


U n -l = 
U 0 ~ 
Wn = 


Vn-t - Vn-t 

h 

U 0 - Uj 

z h 

«„_/ - «„ 


The needed second derivative of the control signals is approximated 

ys - s/2 


ys = 


Vn-S = 


u t = 


Vn — 3 J/n— 4 

h 

yt - i/3 


Un-Z = 
U t = 


2 In — 3 J/n — 2 

h 

Ut - U t 


U n -1 = 


U 0 = 


U n -8 - U n -i 

h 

Uo — Ut 


U n = 


U n -t - U n 


as: 


The Matlab program that solves the task for system 3 using Gaussian elim- 
ination is displayed in mprl51knovel.m. 
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7 The Test 

The test was effected by forcing the current system to run through points 
situated on the curve we wanted to track. This denotes that the trajectory 
has total freedom between the fixed coordinates as long as it passes through 
the test curve dots. How close the system followed the test curve was mea- 
sured at five additional points between each two fixed coordinates. All these 
extra measuring points, along the curve, were added by their absolute value. 
The sum of ail these measurements is called total-error. 

Point-error shows how precisely the system tracks the fixed coordinates. 
During all the trials, this value always been zero, i.e. perfect tracking. The 
only fixed coordinate that the system does not run through is the end point. 
It is due to the lack of constraints there and its divergence is measured by 
End-point-error. 

7.1 Test curves 

Three different kinds of test curves with diverse characteristics were used. 
The two standard curves are one period of the sharply curving sine function 
and the soft curving function, the hyperbolic tangent. A discontinuous step 
function was also used in the test, (see below). 


Or* dim. fraj: mpr12tdar0.m lambda * -1 n « 10 



Figure 1: The step function tracked by system 1 with A=0. 
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One dim. tra|: mpr121d«0jn lambda * -1 n - 26 



Figure 2: The sine function tracked by system 1 with A=0. 


One dim. traj: mpr121def0.m lambda - -1 n - 20 



Figure 3: The tangent hyperbolic function tracked by system 1 with A=0. 
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8 Results 

In this chapter we are going to look at the output from our systems. The 
following figures are exactly like those shown on the computer screen. 

At the top left of the graph the program used is exhibited. At the top right 
are the values of XI , 11, A2, 12, and e, ep, shown together with the number 
of points, n, used to determine the test function. The scale of the y-axis are 
indeterminable but can give a hint about the ratio between variables of the 
same kind. Which quantity the plot gives information about is displayed to 
the left of the axis. At the bottom is always the time axis, scaled in seconds, 
displayed jointly with the calculated errors, see chapter 7. 

At first the difference between system 1 and the other two systems will 
be shown. The two systems can be represented by system 3 since they have 
about the same behavior for this choice of parameters. Notice the sharper 
look of system l’s control signal u and acceleration, the latter has less maxi- 
mal deviation in both graphs. System 1 can not be affected in the same way 

o. 

o. 

0. 

0 . 

ACC 
CSu 

- 0 . 

- 0 . 

- 0 . 


One dim. traj: mpf121dert).m lambda - -1 n * 26 



Figure 4: Sine tracked by system 1. 

as the other two systems. Due to this it is not so very interesting and will 
therefore from now on be omitted. 



8 RESULTS 


41 


On« dim. tr»J: mpr151knov«l.m II - -1 C--1«p-0n-26 



Figure 5: Sine tracked by system 3. 

System 2’s behavior when tracking the soft curving function tangent hy- 
perbolic is shown below. The first graph shows the smooth behavior for the 
system when e = 0. The other two indicate a more uncontrolled fluctuation 
for an e/ 0. The second and the third test are done with different e so the 
received maximum deviation for the acceleration felt by the passengers are 
of the same magnitude. The acceleration is at the bottom except when it is 
oscillating heavily as in the last two plots. 

Opposite the analysis made in chapter 5 it seems as though system 2 is 
not as easy to handle even for an t > 0 and its oscillations for t < 0 are 
apparent. 

In two plots some parts of the curves are omitted. This is to prevent the 
large deviation of the acceleration at the endpoints to suppress other impor- 
tant information. However, this makes a correct error estimation impossible 
and the displayed values on the total error for these plots are wrong. The 
true value on the total error is 0.0997 for figure 7 and 0.1967 for figure 9. 
The last mentioned plot shows the influence of a changed value on A 2. It 
increases the magnitude and shifts the oscillating acceleration to an earlier 
time interval. For both systems it is true that XI affects the behavior much 
more than what A 2 does. 
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One dim. traj: mprl 41 knovel.m II - -0.1 12 » -0.2 ep - 0 n- 20 



Figure 6: Tangent hyperbolic tracked by system 2. 


Ooe dim. traj: mprl 41 knovohm II- -0.1 12 - -0.2 op - -0.01 n - 20 



Figure 7: Tangent hyperbolic tracked by system 2. 
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System 3’s conduct for A 1 = A 2 — -1 and varying values on e, when tracking 
the tangent hyperbolic function, is shown below. 

According to the theoretical discussion in chapter 5, we stated that system 
3 has two zeros in the open right half plane for an e > 0 but only one for an 
e < 0. We see that the system can handle negative e better than positive, 
(see figure 11 and 12). 

The curve that shows the acceleration is mostly beneath the control signal 
u in all graphs. It is also the most oscillating signal in the two last plots. 

The tests using the tangent hyperbolic function are carried out with dif- 
ferent sets of A for each systems. It is therefore not so easy to compare 
the behavior for the actual systems. However, during experiments that are 
not presented in the report, it has been shown that system 3 is more easily 
disturbed for an t ^ 0 than system 2 is. 

Correct total error for figure 11 is 0.2435 and 0.3126 for figure 12. 

Orwdfm. traj: mpn51knoval.m II * -1 12 - -1 op-On-20 



Figure 10: Tangent hyperbolic tracked by system 3. 
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System 3’s characteristics are further examined below when it is applied to 
the sine curve. Figure 5 shows the obtained control signal u and acceleration 
for XI = X2 = -1 and e = 0. The first graph below displays the control signal 
w for the same system and parameters. This is to exhibit the calculated 
signals for a smooth case so we have something to compare with when it is 
getting rough. 


On« dim. tr*J: rrprl 51krwv#Lm 11. -1 t2--1«p-0n-26 



Figure 13: Sin tracked by system 3. 

We receive some other signals in the following figures when system 3 is run by 
an e 7^ 0. Observe the twisted trajectory in figure 15 with an accompanying 
large value on the total error. Notice also the magnitude on the control signal 
w and its third derivative, shown in figures reffi:16 and 17. The system seems 
to be more sensitive for an e > 0 than for an e < 0, this is probably due to 
the circumstance that it has two zeros in the open right half plane in the 
first case but only one zero in the second case. The oscillating acceleration 
is shifted to a later time interval and has there a larger magnitude for a 
negative e. 
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One dm. traj: mpr151knov#J.m II - -1 G - -1 ep - 0.002 n - 28 



Total error : 0.2617 Point error :0 Graph : 8 Time t 


Figure 16: Sine tracked by system 3. 



Total error : 0.2817 Point error : 0 Graph: 11 Timet 


Figure 17: Sine tracked by system 3. 
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On* dfrn. trmj: mpr15Hmcv*lm II - -1 12 - -1 *p - -0.001 n - 28 



Figure 18: Sine tracked by system 3. 


On* dkn. traj: mpr151knov»Lm II - -1 12 - -1 ep - - 0.001 n - 28 



Figure 19: Sine tracked by system 3. 
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Applying the step function to system 3 gives that the very large control 
signals u can be reduced by an e ^ 0 to the cost of a larger toted error. 


On* dUrt traf: mprt 51 knov*Lm N - -0.1 S - -0 l 2 *p - 0 n - 10 



Figure 20: The step function tracked by system 3. 


On* dim. traf mpr151knov*Lm H - -0.1 12 ■ -02 *p« -0.00255 n - 10 



Figure 21: The step function tracked by system 3 
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On# dim. tr«j: mpr151knov«Cm II- -0.1 (2 - -0-2«p - 0.0066 o - 10 



Figure 22: The step function tracked by system 3. 

At last we will take a look at system 2 applied to the sine function with A 1 
— \2 = -1 and for a very particular choice of e. Figures 23 and 24 show 
an oscillating but quite normal behavior for this system. When examining 
figure 24 which is run with a slightly changed e it looks about as the other 
two until the scale on the y-axis is observed. The first two graphs use an t = 
0.004186 respectively e = 0.004188. The value of c for the last three figures 
is 0.00418702471143, which gave the largest oscillatory motions. It seems 
that we have found a set of parameters that brings system 2 in to resonance. 
When changing the value on t we might affect the transfer function so that 
the frequency of the actual in signal w is the peak frequency; see figure 
27. In such a case we will receive an increased amplification of the output 
signal. Notice that this phenomenon only appears when tracking the sine 
test function. 
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On© dim. traj: rTprt41knov©lm II - -1 XZ - -1 ap • 0.004106 n - 26 



Figure 23: Sine tracked by system 2. 


On© dim. traj: rrpii41knov©lnrc II - -1 12 - -1 ep - 0.004188 n - 20 



Figure 24: Sine tracked by system 2. 
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Figure 25: Sine tracked by system 2. 



Figure 26: Sine tracked by system 2. 
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Figure 27: Sine tracked by system 2. 


8.1 Surprising Results 

Looking at figure 4 or 5 one might be surprised at the appearance of the 
control signal u and the acceleration. Even if the sine function is curving it 
should not cause such peaks in the graphs. The reason for this is as follows. 
Consider the system 



From the given system we have that ij = xg and that £g = u. If we use 
splines and force the system to track the function f(t) it implies that the 
first element in the state vector, x;, equals the function, i.e. y(tj = f(t). 
This denotes that x s (t) = f(t) and that the control signal u(t ) = /(<). The 
velocity in figure 28 should by our theory have something in common with 
the derivative of the curve f(t) = sin t. The derivative /(t) = cos t , so in the 
beginning of the first region when f(t) = -1, the graph is shifted upwards 
and scaled, the velocity should be zero. When f(t) becomes zero the velocity 
ought to reach its maximum and then decline. For the other half we have the 
reversed situation and should therefore have an inverted curve. This behavior 
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can be recognized in the first graph and a similar reasoning for the control 
signal u equals f(L) = — sin t gives the calculated control signal in figure 4. 
Here A ^ 0 so the graph is shifted to the right. So the splines do not only 
try to follow the specified trajectory they also approximates its derivatives. 
We can thus effect a perfect trajectory of the test curve but the prize we pay 
is huge values on the control signals and accompanying acceleration. The 
behavior described makes it impossible to realize a smooth aircraft control 
using splines. 

This very simplified heuristic description of the phenomenon will be com- 
pleted in a paper written at Texas Tech University, Lubbock, USA, by Horn 
Professor Clyde Martin, PhD, and Assistant Professor Zhimin Zhang, PhD. 


On# dim. traj: mpr121d#rt).m lambda -On -26 



Figure 28: Sine tracked by system 1. 
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9 Resume in Swedish 

Inledningsvis var varan avsikt att ta fram styrlagar for en flygplans modell 
sa att den flogs sa behagligt for passagerarna som mojligt. Passagerar kora- 
forten var det aJlra mest vasentliga sa vi utvecklade tva kontroll lagax som 
anvande derivatorna av den pilot inducerade styrningen u. Detta ger inte den 
energi snalaste insignalen men tar bort de varsta topparna hos styrsignalen 
och medfoljande acceleration. 

Programvaran som anvants inkluderar Matlab och Maple for berakningar 
och Latex som ordbehandlings program. Andra halften av rapporten bestar 
av Matlab program och som avslutas med en referenslista. 

Vi ansatte den vanliga endimensionella behandlings proceduren och im- 
plementerade systemen i Matlab. Ett huvud program, med tillhorande hjalp 
program, for varje kontroll lag. 

De givna systemen analyserades ur uppnabarhets och stabilitets synpunkt 
vilket resulterade i en bedomning att de var stabila men inte garanterade 
insignal-utsignal stabilitet. Se kapitel 5. 

En “spline” ar en kurva till ett n:te gradens polynom vilken ar forenat med 
liknande polynoms kurvor i respektive andpunkt. I varje forenings punkt har 
funktionerna sina n-1 forsta derivator gemensammma. Detta ger en kurva 
som av ogat tycks vara helt homogen men som i sjalva verket bestar av 
ett antal sammankopplade delar. Kapitel 6 behandlar “splines” och den 
anvandna kontroll teorin. 

Huvud resultatet var att vi inte kunde ta fram mjuka styrlagar eftersom 
“splinsen” inte bara forsoker approximera test kurvan utan aven tar hansyn 
till dess derivator. Genom systemet paverkas aven styrsignalen och vi erhaller 
omojligt stora styrsignaler och accelerationer. Noggranheten i foljningen ar 
alltid exemplarisk. Kapitel 8 behandlar rapportens huvudresultat. 
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10 Mat lab Programs 

function x ■ mprl21der0(spc_plot ,t ,n,cleargr, pointfen .lambda) 

VI. THIS PROGRAM CALCULATES CONTROLLAWS FOR A ONE 
DIMENSIONAL TRAJECTORY VI, 

VI. spc.plot DETERMINES WHICH GRAPH TO BE DISPLAYED, 

CHOSE AN INTEGER =<5 VI. 

VI. t IS THE TIME PERIOD FOR WHICH THE SYSTEM 
IS TO BE CONTROLLED VI. 

VI. n ARE THE NUMBER OF POINTS AT THE SPECIFIED 
TRAJECTORY, CHOSE t/n>l/10 VI. 

VI. cleargr CLEARS THE CURRENT WINDOW, CHOSE 1 OR 0 VI. 

VI. point fen SPECIFIES THE TRAJECTORY VI. 

VI. lambda AFFECTS THE INSTABILITY OF THE SYSTEM VI. 


global A B h t n 
VI. THE SYSTEM VI. 

A=[ 0 1 ; 

0 lambda] ; 

B=[ 0 ; 

1 ]; 

C= [ 1 0 ] ; 

VI. FUNCTION pointfen DETERMINES THE SPECIFIED TRAJECTORY VI. 

VI. NECESSARY FOR THE COMPOUND FUNCTION pointsinl2 VI. 

if pointfcn== , pointsinl2' ; 
t=5 . 2 ; 
n=26 ; 
end; 


R=feval (pointfen) ; 
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XX CALCULATION OF THE INTEGRAL FROM 0 TO h XX 

m=48 ; XX NUMBER OF POINTS BETWEEN INTERPOLATIONS, 

CHOSE A MULTIPLE OF 6 XX 
mp=m/6 ; XX DETERMINES THE PRECISION IN THE 

SPLINE APPROXIMATION XX 

tol=le-08; XX THE NUMERIC ERROR TOLERANCE XX 

Mtau( : , 1 :2)=zeros(2) ; 
tau=0 ; 
for j=l:m 
oldtau=tau; 
tau=oldt au+h/m; 

Mtau( : ,2*j+l :2*j+2)= quad812mod( ' integrand' ,oldtau,tau,tol) 
+ Mtau(: ,2*j-l:2*j) ; 
end; 

M=Mtau( : ,2*m+l :2*m+2) ; 

e_Ah=expm(-A*h) ; 

Minv=inv(M) ; 

ZZ=Minv*e_Ah; 

WW=e_Ah’*ZZ+Minv; 

WL=[WW(2,2)] ; XX PARTITIONING MATRICES XX 
ZLU=[ZZ(2,2)] ; 

ZLL=[ZZ(2,2)] ; 

WR= [WW(2, 1)] ; 

ZRU=[ZZ(2,1)] ; 

ZRL= [ZZ(1 ,2)] ; 

Ulu= [Minv(2,2)] ; 

Uru= [Minv(2 , 1)] ; 

XX FORMING OF THE RIGHT HAND SIDE OF THE BLOCKDI AGONAL 
SYSTEM XX 

for i=2:n 

Omega(i , 1)=-ZRL*R(1 , i-1) +WR*R ( 1 , i ) -ZRU*R ( 1 , i+ 1 ) ; 
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end; 

Omega(l , 1)= Uru*R(l,l) - ZRU*R(1,2); 

Omega(n+l , 1)=ZRL*R(1 ,n) + (Uru-WR)*R(l,n+l) ; 

7.7, FORMING OF THE LEFT HAND SIDE OF THE BLOCKDIAGONAL 
SYSTEM 7 . 7 . 

for i=2:n 

DD(i,i-l)=ZLL; 

DD(i,i)=-WL; 

DD(i,i+l)=ZLU; 

end; 

DD(1 , l)=(lambda-Ulu) ; 

DD(1,2)=ZLU; 

DD(n+l,n)=-ZLL; 

DD(n+l J n+l)=(lambda+WL-Ulu) ; 

DD=sparse(DD) ; •/.*/. SQUEEZING OUT ALL ZERO ELEMENTS 

FROM MATRIX DD 7 . 7 . 

7 . 7 . GAUSSELIMINATION TO PRODUCE AN UPPER 
TRIANGULAR SYSTEM •/.*/. 

for i=l:n 

zd=DD(i+l,i)/DD(i,i) ; 

DD(i+l , i)=DD(i+l , i)-zd*DD(i , i) ; 
DD(i+l,i+l)=DD(i+l,i+l)-zd*DD(i,i+l) ; 

Omega(i+l , l)=0mega(i+l , l)-zd*Omega(i , 1) ; 
end; 

7 . 7 . BACKSUBSTITUTION TO SOLVE FOR THE XVEL */.*/. 

xvel(l ,n+l)=Omega(n+l , l)/DD(i+l , i+1) ; 
for i=n:-l:l 

xvel(l , i) = (Omega(i , l)-DD(i,i + l)*x ve l(l»i. + l))/DD(i, i) ; 
end; 
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’/.•/. MAKING OF THE STATE VECTOR '/.*/. 
for i=0:n 

x(: ,i+l)=[[R(l,i+l)] 

[xvel (1 , i+l)]] ; 

end; 

if cleargr 
elf ; 

hold on; 
grid on; 
end; 


*/.*/. PLOTTING OF THE CALC/SPEC TRAJECTORY, VELOCITY, 
CONTROLS IGNAL AND ACCELERATION 1 . 1 . 

plotadm=0 ; 
while spc_plot 
if plotadm 

if spc_plot<6 
figure; 
hold on; 
grid on; 

title(['0ne dim. traj : mprl21der0.m lambda = ' , 
num2str (lambda) , ’ n = ' ,num2str(n)] ) 
xlabel( ['Total error : ' ,num2str(Total_error) , ' 
Pointerror : ’ ,num2str(Point_error) , * Graph : 

’ ,num2str(spc_plot) , ' Time t']) 
for k=0:n 

plot (k*h,x(l ,k+l) , 'o’) 
end; 
end; 
end; 

breakadm=0; 

f adm=0; 1 . 1 . NORMALLY fadm=0, NECESSARY FOR WRITING OF 
TEXTS IN THE PLOT 1 . 1 . 
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for j=0:m 

eAtau=expm(A*j*h/m) ; 

7 , 7 . CHOOSE e.g 2 :n-3 TO AVOID PROBLEMS AT THE ENDPOINTS 7 . 7 . 
for i»fadm:n-l 

entry( : , i+l)=eAtau*(x( : , i+l)+Mtau( : ,2*j+l :2*j+2)* 
Minv*(e_Ali*x( : ,i+2)-x(: ,i+l))) ; 
csignvec ( : , i+ 1 ) =B ’ *expm(-A ‘ * j *h/ m) *Minv* 

(e_Ah*x( : , i+2) -x( : , i+1) ) ; 
if j==0 

if i==f adm; 

entryl=entry (1 ,f adm+1) ; 
entry2=entry(2,fadm+l) ; 
csignvecl=csignvec(l ,fadm+l) ; 
end; 
end; 

if rem(j,mp)==0 
if j<=m-mp 

traj ect(l , j/mp*n+(i+l))=entry(l,i+l) ; 
end; 
end; 
if j==m 


if i==n-l 

traject(l,6*n+l)=entry(l,i+l) ; 
end; 
end; 

if spc_plot==l 

plot (i*h+j *h/m, entry (1 , i+1) , ' . ') 
plot(i*h+j*h/m,csignvec(l,i+l) , ’ . 

7 . 7 . ACCELERATION 7 . 7 . 

plot (i*h+j*h/m,lambda*entry(2, i+l)+ 
csignvec(l , i+1) , ' . ’ ) 
elseif spc_plot==2 

plot(i*h+j*h/m,entry(l,i+l) , ’ . ') 
plot ( i*h+ j*h/m, entry (2, i+1) , ' . *) 

plot(i*h+j*h/m,csignvec(l,i + l) , ’ • ’) 

7 . 7 . ACCELERATION •/.*/. 

plot (i*h+ j *h/m , lambda*entry (2 , i+1) + 


7 . 7 . TRAJECTORY 7 . 7 . 
7 . 7 . CONTROL u 7 . 7 . 


7 . 7 . TRAJECTORY 7 . 7 . 
7 . 7 . VELOCITY 7 . 7 . 

7 . 7 . CONTROL u 7 . 7 . 
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csignvec(l , i+1) , ’ . ’ ) 
elseif spc_plot==3 

plot(i*h+j*h/m,entry(l,i+l) , ’ . ’) 

XX ACCELERATION */.*/. 

plot (i*h+j*h/m, lambda*entry(2,i+l)+ 
csignvec(l , i+1) , ' . ') 
elseif spc_plot==4 

plot (i*h+j*h/m, entry (1 , i+1) , ’ . ‘ ) 
plot (i*h+j*h/m, entry (2, i+1) , ' . ' ) 
elseif spc_plot==5 

plot(i*h+j*h/m,entry(l,i+l) , ’ . ') 
plot(i*h+j*h/m,csignvec(l,i+l) , ’ . *) 
else 

disp(' ') 

disp( ’ NOT A VALID CHOICE ’) 
breakadm=l ; 
end; 

if breakadm 
break; 
end; 
end; 

if breakadm 
break; 
end; 
end; 

if spc_plot<6 
ax = ax is ; 
ax2=ax (1,2) ; 
if ax2<1.5 
xpos=3/4*0.2; 
elseif ax2>=5 
xpos=3/4; 
else 

xpos=0.25; 

end; 

if spc_plot<3 
ax4=ax(l ,4) I 


VI. TRAJECTORY VI. 


VI. TRAJECTORY VI. 
U VELOCITY VI. 

VI. TRAJECTORY VI. 

y.y, control u y.y. 
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if ax4<l 
ax4=ax4*10 ; 
if ax4<l 

ax4=ax4*10 ; 
end; 
end; 

if ax4>10 
ax4=ax4/10; 
if ax4>10 
ax4=ax4/10 ; 
end; 
end; 

if rem(ax4,2)==0 
yadd-(ax(l ,4)/4) *1/5 ; 
else 

yadd=(ax(l ,4)/3)*l/5 ; 
end; 

ypos=lambda*entry2+csignvecl ; 
text(-xpos ,ypos, ['ACC’]) ; 
if abs(ypos-csignvecl)>yadd 

text(-xpos,csignvecl, ['CS u J ]); 
elseif ypos-csignvecl>0 

text (-xpos,csignvecl-yadd, ['CS u']) ; 
else 

text (-xpos,csignvecl+yadd, ['CS u’]) ; 
end ; 

if spc_plot==2 

text(-xpos,entry2, ['VEL'] ) ; 
end; 
end; 

if spc_plot==3 

ypos=lambda*entry2+csignvecl ; 
text (-xpos ,ypos , ['ACC']) ; 
end; 

if spc_plot==4 

text (-xpos ,entry2 , ['VEL'] ) ; 
end; 
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if spc_plot==5 

text(-xpos,csignvecl, [’CS u']); 
end; 
end; 


'/.*/. ESTIMATION OF THE ACCURACY IN THE SPLINE 
APPROXIMATION 7.’/. 


if plotadm==0 

[Total_error , Average_error , Point_error,End_point_error] = 
spl ine_error (po intf cn , R , traj ect ) ; 

Total_error 

Average_error 

Point_error 

End_point_error 

title(['0ne dim. traj: mprl21der0.m lambda * ’ ,num2str 
(lambda) , ' n = ’ ,num2str(n)] ) 

xlabel( ['Total error : ’ ,num2str(Total_error) , ' Point 
error : ' ,num2str(Point_error) , ' Graph : 
num2str(spc_plot) , ' Time t']) 
if spc_plot<6 
for k=0:n 

plot (k*h,x(l ,k+l) , ' o ' ) 
end; 
end; 
end; 

plotadm=plotadm+l ; 
disp(' ') 

disp( ' WOULD YOU LIKE TO SEE ANOTHER GRAPH OF THE CURRENT 
SYSTEM AND ITS TRAJECTORY ? ') 
disp(' ') 

disp( ' FINISH THE PROGRAM : 0 ') 

disp( ' DISPLAY THE TRAJECTORY, CONTROL u AND 

ACCELERATION : 1 ') 

disp( ' DISPLAY THE TRAJECTORY, VELOCITY, CONTROL u AND 
ACCELERATION : 2 ') 
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disp( ' DISPLAY THE TRAJECTORY AND ACCELERATION : 3 ’) 
disp( J DISPLAY THE TRAJECTORY AND VELOCITY : 4 ') 
disp( ' DISPLAY THE TRAJECTORY AND CONTROL u : 5 ') 
spc_plot = input (' MAKE YOUR CHOICE : ’); 
end; 


clear global; 
for k=2:plotadm 
delete(k) ; 
end; 
end; 


function [Q,cnt] = quad812mod(funfcn,a,b,tol) 

^Alteration of the original matlab toolbox program. 

XQUAD8 Numerical evaluation of am integral, higher order 
*/, method. Q = QUAD8(’F’ ,A,B,T0L) approximates the 
7 , integral of F(X) from to B to within a relative error 

7, of TOL. 'F' is a string containing the name of the 

7, function. The function must return a 2*2-matrix 
7 , output value if given an input value. 

*/, Q = Inf is returned if am excessive recursion level 
7 , is reached indicating a possibly singular integral. 

7 , QUAD8 uses an adaptive recursive Newton Cotes 8 pamel 

7, rule . 

7, Cleve Moler, 5-08-88. 

7, Copyright (c) 1984-94 by The MathWorks, Inc. 

7% [Q.cnt] = quad8(F,a,b,tol) also returns a function 

7 , evaluation count. 

7, Top level initialization, Newton-Cotes weights 

w = [3956 23552 -3712 41984 -18160 41984 -3712 23552 
3956] /14175 ; 

x = a + (0:8)*(b-a)/8; 


•/. set up function call 
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for i=x 

y = [y feval(funfcn.i)] ; 
end; 

X Adaptive, recursive Newton-Cotes 8 panel quadrature 
QO = zeros(2) ; 

[Q,cnt] = quad812stpmod(funfcn,a,b,tol,0,w,x,y,Q0) ; 

cnt = cnt + 9; 

end; 


function [Q,cnt] = quad812stpmod(FunFcn,a,b,tol,lev, 

w,x0,f0,Q0) 

XAlteration of the original matlab toolbox program. 
XQUAD8STP Recursive function used by QUAD8. 

X [Q,cnt] * quad8stp(F,a,b,tol,lev,w,f ,Q0) tries to 
X approximate the integral of f(x) from a to b to 
X within a relative error of tol. F is a string 
X containing the name of f . The remaining arguments 
X are generated by quad8mod or by the recursion. 

X lev is the recursion level. 

X w is the weights in the 8 panel Newton Cotes formula. 

X xO is a vector of 9 equally spaced abscissa is the 

X interval. 

X fO is a matrix of the 9 function values at x. 

X QO is an approximate value of the integral. 

X Cleve Moler, 5-08-88. 

X Copyright (c) 1984-94 by The Mathtforks, Inc. 

LEVMAX = 10; 

X Evaluate function at midpoints of left and 
right half intervals, 
x = zeros(l , 17) ; 
x(l :2: 17) = xO; 

x(2 :2 : 16) = (x0(l:8) + x0(2:9))/2; 
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f ( : , 1 : 2)= fO(: ,1:2) ; 

for i=l : 8 

f ( : ,4*i-l :4*i) = feval(FunFcn,x(2*i)) ; 
f (: ,4*i+l:4*i+2) = (f 0( : ,2*i+l : 2*i+2) ) ; 
end; 

X Integrate over half intervals, 
h = (b-a)/16; 

Q1=0 ; Q2=0 ; 

for i=l : 9 

Q1 = Q1 + h*w(i)*f ( : ,2*i-l :2*i) ; 
q2 = Q2 + h*w(10-i)*f ( : ,35-i*2:36-i*2) ; 
end; 

q = qi + q2; 

X Recursively refine approximations, 
if norm(q - qo) > tol*norm(q) & lev <= LEVMAX 
c = (a+b)/2; 

[qi.cntl] = quad812stpmod(FunFcn,a,c,tol/2,lev+l, 

w,x(l:9),f(:,l:18),qi); 
Cq2,cnt2] * quad812stpmod(FunFcn,c,b,tol/2,lev+l, 

w,x(9:17),f(:, 17:34), q2); 

q = qi + q2; 

cnt = cnt + cntl + cnt2; 

end 

end; 
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function x = mprl41knovel(spc_plot,t,n,cleargr, pointfcn, 
lambdal , lambda2 , ep) 

7 . 7 . THIS PROGRAM CALCULATES CONTROLLAWS FOR A ONE 
DIMENSIONAL TRAJECTORY 7.7. 

7.7. spc.plot DETERMINES WHICH GRAPH TO BE DISPLAYED, 
CHOSE AN INTEGER =<9 7.7. 

7.7. t IS THE TIME PERIOD FOR WHICH THE SYSTEM IS 
TO BE CONTROLLED 7.7. 

7.7. n ARE THE NUMBER OF POINTS AT THE SPECIFIED 
TRAJECTORY, CHOSE t/n>l/10 7.7. 

7.7. cleargr CLEARS THE CURRENT WINDOW, CHOSE 1 OR 0 7.7. 

7.7. pointfcn SPECIFIES THE TRAJECTORY 7.7. 

7.7. lambdal AFFECTS THE INSTABILITY OF THE SYSTEM 7.7. 

7.7. Iambda2 ALSO EFFECTS THE STABILITY OF THE SYSTEM, 
CHOSE 11~=12 7.7. 

7.7. ep<>0 PUTS A ZERO IN THE TRANSFERFUNCTION 7.7. 


global A B h t n 

7.7. THE SYSTEM 7.7. 

A=[ 0 1 0 0; 

0 lambdal 1 0; 

0 0 0 1 ; 

000 lambda2] ; 

B=[ 0 ep 0 1] ’ ; 

C= [ 1 0 0 0] ; 

7.7. FUNCTION pointfcn DETERMINES THE SPECIFIED TRAJECTORY 7.7. 

7.7. NECESSARY FOR THE COMPOUND FUNCTION pointsinl2 7.7. 

if pointfcn" ’pointsinl2' ; 
t=5 . 2 ; 
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n=26; 

end; 


R=feval(pointfcn) ; 

*/,•/. CALCULATION OF THE INTEGRAL FROM 0 TO h VI. 

m=48; •/.*/. NUMBER OF POINTS BETWEEN INTERPOLATION, 

CHOSE A MULTIPLE OF 6 VI. 

mp=m/6; V/. NEEDED FOR fen spline_error THAT 

DETERMINES THE PRECISION IN THE 
SPLINE APPROXIMATION V/. 

tol=le-08; VI. THE NUMERIC ERROR TOLERANCE VI. 

Mtau ( : ,1:4) =zeros (4) ; 
tau=0 ; 
for j = l:m 
oldtau=tau; 
tau=oldtau+h/m; 

Mtau( : ,4*j+l :4*j+4)= quad814mod( ’ integrand’ ,oldtau,tau,tol) 
+ Mtau( : ,4*j-3:4*j) ; 
end; 

M=Mtau( : ,4*m+l :4*m+4) ; 


*/.*/. FORMING OF THE MATRICES FOR THE BLOCKDIAGONAL SYSTEM */.•/, 

e_Ah=expm(-A*h) ; 

Minv=inv(M) ; 

ZZ=Minv*e_Ah; 

WW=e_Ah ' *ZZ+Minv ; 

•/.*/. CONTINUOUS CONTROLLAW 

WLDO= [ep*WW (2 , 2) +WW (4 , 2) ep*WW(2 ,3)+WW(4,3) ep*WW(2,4)+ 
WWC4.4)]; 

ZLUDO= [ep*ZZ (2 , 2) +ZZ (4 , 2) ep*ZZ(2,3)+ZZ(4,3) ep*ZZ(2,4)+ 
ZZ(4,4)] ; 
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ZLLD0= [ep*ZZ(2 ,2) +ZZ(2 ,4) ep*ZZ(3,2)+ZZ(3,4) ep*ZZ(4,2)+ 
ZZ(4,4)3 ; 

WRD0=[ep*WW(2,l)+WW(4,l)] ; 

ZRUD0= [ep*ZZ(2 , 1)+ZZ(4, 1)] ; 

ZRLD0= [ep*ZZ ( 1 , 2 ) +ZZ ( 1 , 4) ] ; 

7 . 7 . CONTINUOUS FIRST DIFFERENTIAL OF CONTROLLAW 7 . 7 . 

WLD1= [ep*WW (1,2) +ep*lambdal*WW (2 , 2) +WW (3 , 2) +lambda2*WW (4,2) 
ep*WW (1,3) +ep*lambdal*WW (2 , 3) +WW (3 , 3)+lambda2*WW (4,3) 
ep*WW (1,4) +ep*lambdal *WW (2 , 4) +WW (3,4) +lambda2*WW (4,4)] ; 

ZLUD1= [ep*ZZ(l , 2)+ep*lambdal*ZZ (2 , 2)+ZZ(3 , 2)+lambda2*ZZ(4, 2) 
ep*ZZ (1,3) +ep*lambdal*ZZ (2,3) +ZZ (3 , 3) +lambda2*ZZ (4 , 3) 
ep*ZZ (1,4) +ep*lambdal *ZZ (2 , 4) +ZZ (3 , 4) +lambda2*ZZ (4,4)] ; 

ZLLD1= [ep*ZZ(2 , l)+ep*lambdal*ZZ(2 ,2)+ZZ(2 ,3)+lambda2*ZZ(2 ,4) 
ep*ZZ (3,1) +ep*lambdal*ZZ (3 , 2) +ZZ (3 , 3) +lambda2*ZZ (3 , 4) 
ep*ZZ (4,1) +ep*lambdal *ZZ (4 , 2) +ZZ (4 , 3) +lambda2*ZZ (4,4)] ; 

WRD1= [ep*WW(l , l)+ep*lambdal*WW(2, 1)+WW(3, l)+lambda2*WW(4, 1)] ; 

ZRUD1= [ep*ZZ(l , 1) +ep*lambdal*ZZ(2 , 1)+ZZ(3 , l)+lambda2*ZZ(4, 1)] ; 

ZRLD1= [ep*ZZ (1,1) +ep*lambdal*ZZ ( 1 , 2) +ZZ (1,3) +lambda2*ZZ (1,4)]; 

7 . 7 . CONTINUOUS SECOND DIFFERENTIAL OF CONTROLLAW 7 . 7 . 

WLD2= [ep*lambdal*WW (1,2)+ (ep*lambdal ‘2+1) *WW(2 , 2) +lambda2* 
WW(3, 2)+lambda2‘2*WW(4,2) ep*lambdal*WW(l ,3)+ 
(ep*lambdal‘2+l)*WW(2 ,3)+lambda2*WW(3,3)+lambda2“2* 

WW (4 , 3) ep*lambdal*WW (1,4)+ (ep*lambdal ‘2+1) *WW (2 , 4) + 
lambda2*WW (3 , 4) +lambda2‘ 2*WW (4,4)] ; 

ZLUD2= [ep*lambdal*ZZ(l , 2)+(ep*lambdal‘2+l)*ZZ(2 ,2)+lambda2* 
ZZ(3,2)+lambda2‘2*ZZ(4,2) ep*lambdal*ZZ(l,3)+ 
(ep*lambdal‘2+l)*ZZ(2 ,3)+lambda2*ZZ(3,3)+ 
lambda2‘2*ZZ(4,3) ep*lambdal*ZZ(l ,4)+(ep*lambdal‘2+l) * 
ZZ(2,4)+lambda2*ZZ(3,4)+lambda2‘2*ZZ(4,4)] ; 

ZLLD2=[ep*lambdal*ZZ(2,l)+(ep*lambdal‘2+l)*ZZ(2,2)+lambda2* 
ZZ(2,3)+lambda2‘2*ZZ(2,4) ep*lambdal*ZZ(3, 1)+ 
(ep*lambdal‘2+l)*ZZ(3,2)+lambda2*ZZ(3,3)+lambda2‘2* 
ZZ(3,4) ep*lambdal*ZZ(4, l)+(ep*lambdal‘2+l)*ZZ(4,2)+ 
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lambda2*ZZ(4,3)+lambda2‘2*ZZ(4,4)] ; 

HRD2=[ep*lambdal*WW(l,l)+(ep*lambdal“2+l)*HW(2,l)+lambda2* 

WW(3, l)+lambda2*2*WW(4, 1)] ; 

ZRUD2=[ep*lambdal*ZZ(l,l)+(ep*lambdal~2+l)*ZZ(2,l)+lambda2* 
ZZ(3,l)+lambda2“2*ZZ(4,l)] ; 

ZRLD2=[ep*lambdal*ZZ(l,l)+(ep*lambdal~2+l)*ZZ(l,2)+lambda2* 
ZZ(l,3)+lambda2~2*ZZ(l,4)] ; 

7.V. DIFFERENTIAL APPROXIMATIONS FOR THE BOUNDARY 
CONDITIONS U 

ydlO»(R(l,2)-R(l,l))/h; 
ydll=(R(l,3)-R(l,2))/h; 
ydl2=(R(l ,4)-R(l ,3) )/h; 
ydln=(R(l ,n+l)-R(l,n))/h; 
ydln_l=(R(l ,n)-R(l,n-l))/h; 
ydln_2=(R(l ,n-l)-R(l ,n-2))/h; 
uO=(ydlO-ydll)/h; 
ul=(ydl l-ydl2) /h ; 
un=(ydln_l-ydln) /h; 
un_l= (ydln_2-ydln_l) /h ; 
udlO=(uO-ul)/h; 
udln= (un_l-un) /h ; 

XO= [yd 10 ; uO; udlO] ; 7.7, 3/4 of the the first state vector 7.7. 
Xn= [ydln; un; udln] ; */.•/, 3/4 of the the last state vector 7.7. 

7 . 7 . FORMING OF THE RIGHT HAND SIDE OF THE 
BLOCKDIAGONAL SYSTEM */.*/. 

j*i; 

for i=2:n 

Omega (j , 1)=-ZRLD0*R(1 , i-l)+WRD0*R(l , i)-ZRUD0*R(l , i+1) 

Omega (j , 1)=-ZRLD1*R(1 , i-l)+WRDl*R(l , i) -ZRUD1*R(1 , i+l) 

Omega (j , l)=-ZRLD2*R(l , i-l)+WRD2*R(l , i)-ZRUD2*R(l,i+l) ;j=j+l; 
end; 

Omega (1 , l)=Omega(l , l)-ZLLDO*XO ; 

0mega(2, l)=0mega(2, 1)-ZLLD1*X0; 
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Omega (3 , l)=0mega(3, 1)-ZLLD2*X0 ; 
0mega(3*n-5, l)=0mega(3*n-5, l)-ZLUDO*Xn; 
Omega(3*n-4, l)=0mega(3*n-4, l)-ZLUDl*Xn; 
0mega(3*n-3, l)=0mega(3*n-3, l)-ZLUD2*Xn; 

'/,*/, FORMING OF THE LEFT HAND SIDE OF THE 
BLOCKD I AGONAL SYSTEM %X 

for i=l:3 

DD(1 , i)=-WLDO(l ,i) ; 

DD(2,i)— WLDl(l.i) ; 

DD(3, i)=-WLD2(l ,i) ; 
if n>2 

DD(1 , i+3)=ZLUD0(l , i) ; 

DD(2 , i+3)=ZLUDl (1 , i) ; 

DD(3 , i+3)=ZLUD2(l , i) ; 
end; 
end; 

for i=2:n-2 
for j =— 1 : 1 

DD (3*i-2 , 3*i-4+j ) =ZLLD0 ( 1 , j +2) ; 

DD (3*i-2 ,3*i-l+j )=-WLD0 (1 , j +2) ; 
DD(3*i-2,3*i+2+j)=ZLUD0(l, j+2); 
DD(3*i-l,3*i-4+j)=ZLLDl(l, j+2) ; 
DD(3*i-l ,3*i-l+j )=-WLDl (1 , j+2) ; 
DD(3*i-l,3*i+2+j)=ZLUDl(l, j+2); 
DD(3*i ,3*i-4+j)=ZLLD2(l, j+2) ; 
DD(3*i ,3*i-l+j)=-WLD2(l, j+2) ; 
DD(3*i , 3*i+2+j)=ZLUD2(l , j+2) ; 
end; 
end; 
if n>2 

for i=l:3 

DD (3*n-5 , 3*n-9+i) =ZLLD0 ( 1 , i) ; 
DD(3*n-4,3*n-9+i)=ZLLDl (1 , i) ; 

DD (3*n-3 , 3*n-9+i ) =ZLLD2 ( 1 , i ) ; 

DD (3*n-5 ,3*n-6+i) =-WLD0 C 1 , i) ; 
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DD(3*n-4,3*n-6+i)=-WLDl(l,i) ; 
DD(3*n-3,3*n-6+i)=-WLD2(i,i) ; 
end; 

end; 

DD=sparse(DD) ; 7.7. SQUEEZING OUT ALL ZERO ELEMENTS 

FROM MATRIX DD 7.7. 

7.7. GAUSSELIMINATION TO PRODUCE AN UPPER TRIANGULAR 
SYSTEM 7.7. 

for k=2:3:3*n-7 
for l=k:k+4 

if DD(k-l ,k-l) ~=0 

zd=DD(l,k-l)/DD(k-l ,k-l) ; 

DD(l , : )=DD(1 , : )-zd*DD(k-l , : ) ; 

Omega(l, l)=Omega(l, l)-zd*Omega(k-l , 1) ; 
end; 
end; 

for l=k+l:k+4 
if DD(k,k) "=0 

zd=DD (l,k)/DD(k,k); 

DD(1, :)=DD(1, :)-zd*DD(k, :); 

Omega (1,1) =0mega (1,1) -zd*0mega (k , 1 ) ; 
end; 
end ; 

for l=k+2:k+4 

if DD(k+l,k+l)'»0 

zd=DD(l ,k+l)/DD(k+l ,k+l) ; 

DD(1, :)=DD(1, :)-zd*DD(k+l, :) ; 

Omega(l , 1) =0mega(l , 1) -zd*Omega(k+l , 1) ; 
end; 
end ; 

end ; 

k=3*n-5 ; 

if DD(k,k)'=0 

zd=DD(k+l,k)/DD(k,k) ; 
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DD(k+l, :)=DD(k+l, :)-zd*DD(k, :); 

Omega(k+l , l)=Omega(k+l , l)-zd*Omega(k, 1) ; 
zd=DD (k+2 , k) /DD (k , k) ; 

DD(k+2,:)=DD(k+2, :)-zd*DD(k, :); 

Omega (k+2 , 1) =0mega(k+2 , 1) -zd*Omega(k, 1) ; 
end; 

if DD(k+l,k+l) - =0 

zd=DD(k+2,k+l)/DD(k+l ,k+l) ; 

DD (k+2 , : ) =DD (k+2 , : ) -zd*DD (k+ 1 , : ) ; 
0mega(k+2 , 1) =0mega(k+2 , 1) -zd*Omega(k+l , 1) ; 
end; 


17 , BACKSUBSTITUTION TO SOLVE FOR THE STATEVECTORS '/.•/. 


udl(n-l)=0mega(3*(n-l) , l)/DD(3*(n-l) ,3*(n-l)) ; 
u(n-l)=(Omega(3*n-4, l)-DD(3*n-4,3*(n-l) )*udl(n-l) )/ 
DD(3*n-4,3*n-4) ; 

ydl (n-l)=(0mega(3*n-5 , l)-DD(3*n-5 ,3*(n-l) )*udl(n-l)- 
DD(3*n-5 , 3*n-4) *u(n-l) ) /DD(3*n-5 ,3*n-5) ; 
if n>2 

for k=n-2:-l:l 

udl (k) = (Omega (3*k, l)-DD(3*k,3*k+3)*udl(k+l)- 
DD(3*k,3*k+2)* u(k+l)-DD(3*k,3*k+l)*ydl(k+l))/ 
DD(3*k,3*k) ; 

u(k)=(0mega(3*k-l , l)-DD(3*k-l ,3*k+3)*udl(k+l)- 
DD (3*k-l ,3*k+2) *u(k+l) -DD(3*k-l , 3*k+l) *ydl (k+1)- 
DD (3*k- 1 , 3*k) *ud 1 (k) ) /DD (3*k- 1 , 3*k-l) ; 
ydl (k) = (Omega (3*k-2 , 1 ) -DD(3*k-2 , 3*k+3) *udl (k+1)- 
DD(3*k-2 , 3*k+2) *u(k+l) -DD (3*k-2 , 3*k+l) *ydl (k+1) 
-DD (3*k-2 , 3*k) *udl (k) -DD (3*k-2 , 3*k- 1) *u(k) )/ 

DD (3*k-2 , 3*k-2) ; 
end; 
end; 


V/. MAKING OF THE STATEVECTORS */.% 


x(: f l)-[R(l f l); XO]; 
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x( : ,n+l)= [R(l ,n+l) ; Xn] ; 
for i=l : n-1 

x( : , i+l)=[R(l, i+1) ; ydl(i); u(i); udl(i)]; 
end; 

if cleargr 
elf ; 

hold on; 
grid on; 
end; 


n PLOTTING OF THE CALC/SPEC TRAJECTORY, VELOCITY, 
CONTROLSIGNAL AND ACCELERATION '/.*/, 

plotadm=0 ; 
while spc.plot 
if plotadm 

if spc_plot<10 
figure; 
hold on; 
grid on; 

title(['0ne dim. traj : mprl41knovel .m 11 = 

> ,num2str(lambdal) , ’ 12 = ’ ,num2str(lambda2) , ’ ep « 

' ,num2str(ep) , ’ n = ’ ,num2str(n)] ) 

xlabel([ 'Total error : ' ,num2str(Total_error) , 

' Point error : ' ,num2str(Point_error) , ’ Graph : 
' ,num2str(spc_plot) , ’ Timet']) 
if spc_plot<6 
for k=0:n 

plot(k*h,x(l ,k+l) , 'o') 
end; 
end; 
end; 
end; 

breakadm=0 ; 

f adm=0 ; VI. NORMALLY fadm=0, NECESSARY FOR WRITING 
OF TEXTS IN THE PLOT */.*/. 
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for j=0:m 

eAtau=expm(A*j*h/m) ; 

for i=f adm:n-l •/.*/. CHOOSE e.g 2:n-3 TO AVOID 

PROBLEMS AT THE ENDPOINTS */.*/. 
entry ( : , i+l)=eAtau* (x( : , i+l)+Mtau( : ,4*j+l :4*j+4) * 
Minv*(e_Ah*x( : , i+2)-x( : , i+l) ) ) ; 

csignvec ( : , i+ 1 ) =B ' *expm( -A ' * j *h/m) *Minv* (e_Ah*x ( : , i+2) - 
x ( : , i+l) ) ; 
if j==0 

if i==f adm; 

entryl=entry(l,fadm+l) ; 
entry2=entry (2,fadm+l) ; 
entry3=entry (3,f adm+1) ; 
entry4=entry (4,f adm+1) ; 
csignvecl=csignvec(l ,fadm+l) ; 
end; 
end; 

if rem(j ,mp)==0 
if j <=m-mp 

traject(l , j/mp*n+(i+l))=entry(l,i+l) ; 
end; 
end; 
if j==m 
if i==n-l 

traject(l,6*n+l)=entry(l,i+l) ; 
end; 
end; 

if spc_plot==l 

plot (i*h+j*h/m, entry (1, i+l) , ' . ’) 7.7. TRAJECTORY 7.7. 

plot(i*h+j*h/m,entry(3,i+l) , ’ . ') 7.7. CONTROL u 7.7. 

7.7. ACCELERATION 7.7. 

plot (i*h+ j *h/m , lambdal*entry (2 , i+l) + 
entry(3, i+l)+ep*csignvec(l , i+l) A) 
elseif spc_plot==2 

plot(i*h+j*h/m,entry(l ,i+l) , ’ . ') 7*1. TRAJECTORY 7.7. 

plot (i*h+j*h/m, entry (2, i+l) , ' . ’) 7.*/. VELOCITY 7.7. 

plot (i*h+j*h/m, entry (3, i+l) , ' . ’) 7.7. CONTROL u 7.7. 
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plot (i*h+j*h/m, lambdal*entry(2,i+l)+entry(3, i+l)+ 
ep*csignvec(l , i+1) , ’ . ’ ) 7.7. ACCELERATION 7.7. 

elseif spc_plot==3 

plot(i*h+j*h/m,entry(l,i+l) , ' . ’) 7.7, TRAJECTORY 7,7. 

plot (i*h+j*h/m,lambdal*entry(2,i+l)+entry(3, i+l)+ 
ep*csignvec(l , i+1) 7.7. ACCELERATION VI. 

elseif spc_plot==4 

plot(i*h+j*h/m,entry(l,i+l) , ' . ') VI. TRAJECTORY VI. 
plot (i*h+j*h/m,entry(2, i+1) ') VI. VELOCITY VI. 
elseif spc_plot==5 

plot (i*h+j*h/m , entry (1 , i+1) ,’.' ) VI. TRAJECTORY VI. 
plot(i*h+j*h/m,entry(3,i+l) , ' . ') VI. CONTROL u VI. 
elseif spc_plot==6 

plot(i*h+j*h/m,entry(4,i+l) , ' . ') VI. CONTROL DER 

u-dot VI. 

elseif spc_plot==7 

csignvec( : ,i+l)=B , *expm(-A'*j*h/m)*Minv* 

(e_Ah*x( : , i+2) -x( : , i+1) ) ; 

7.7. CONTROLS IGNAL w •/,*/. 

plot (i*h+j*h/m,csignvec(l , i+1) , ’ . ’ ) 

elseif spc_plot==8 

csdlvec( : J i+l)=B , *A , *expm(-A’*j*h/m)*Minv* 

(e_Ah*x( : ,i+2)-x( : ,i+l)) ; 

7.7. CONTROL DER w-dot 7.7. 
plot(i*h+j*h/m,csdlvec(l,i+l) , 1 . ’) 
if j==0 

csdvecl=csdlvec(l , 1) ; 
end; 

elseif spc_plot==9 

csd2vec( : , i+1) =B' *A ' *A ' *expm(-A ’ * j *h/m) * 
Minv*(e_Ah*x( : , i+2)-x( : , i+1) ) ; 

7.7. CONTROL DERx2 w-dot-dot */.•/. 
plot(i*h+j*h/m,csd2vec(l,i+l) , ' . ') 
if j==0 

csdvec2=csd2vec(l , 1) ; 
end; 
else 
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disp(' ') 

disp( ' NOT A VALID CHOICE ’) 
breakadm=l ; 
end; 

if breakadm 
break; 
end; 
end; 

if breakadm 
break; 
end; 
end; 

if spc_plot<10 
ax=axis ; 
ax2=ax(l ,2) ; 
if ax2<1.5 
xpos=3/4*0.2; 
elseif ax2>=5 
xpos=3/4; 
else 

xpos=0.25; 

end; 

if spc_plot<3 
ax4=ax(l ,4) ; 
if ax4<l 
ax4=ajc4*10; 
if ajc4<l 
ax4=ajc4*10; 
end; 
end; 

if ax4>10 
ax4=ax4/10 ; 
if ajc4>10 
ajc4=ax4/ 10 ; 
end; 
end; 

if rem(ax4 , 2) ==0 
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yadd=(ax(l ,4)/4)*l/5; 
else 

yadd=(ax(l,4)/3)*l/5; 

end; 

ypos=lainbdal*entry2+entry3+ep*csignvecl; 
t ext ( -xpos, ypos, ['ACC'] ) ; 
if abs(ypos-entry3)>yadd 

text(-xpos,entry3, [’CS u']); 
elseif ypos-entry3>0 

text (-xpos, entry3-yadd, ['CS u']) ; 
else 

text (-xpos, entry3+yadd, ['CS u'] ) ; 
end; 

if spc_plot==2 

text(-xpos,entry2, ['VEL']) ; 
end; 
end; 

if spc_plot==3 

ypos=lajnbdal*entry2+entry3+ep*csignvecl; 
text (-xpos ,ypos, ['ACC'] ) ; 
end; 

if spc_plot==4 

text (-xpos, entry2, ['VEL']); 
end; 

if spc_plot==5 

text (-xpos, entry3 , ['CS u']); 
end; 

if spc_plot==6 

text (-xpos, entry4, ['udl']) ; 
end; 

if spc_plot==7 

text (-xpos , csignvecl , ['CS w']); 
end; 

if spc_plot==8 

t ext ( -xpos, csdvecl, ['wdl']) ; 
end; 

if spc_plot==9 
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text(-xpos,csdvec2, ['vd2']) ; 
end; 
end; 


7 . 7 . ESTIMATION OF THE ACCURACY IN THE 
SPLINE APPROXIMATION 7 . 7 . 


if plotadm==0 

[Tot al_error , Average_error , Point _error , 

End_po int _error] =spline_error (po intf cn ,R , traj ect ) ; 

Total_error 

Average_error 

Point_error 

End_point_error 

title(['0ne dim. traj: mprl41knovel .m 11 — 

' ,num2str(lambdal) , ’ 12 = ’ ,num2str(lambda2) , 

’ ep = ’ ,num2str(ep) , ' n * ’ ,num2str(n)]) 
xlabel( ['Total error : ' ,num2str(Total_error) , 

' Point error : ’ ,num2str(Point_error) , ' Graph : 
' ,num2str(spc_plot) , ' Time t']) 
if spc_plot<6 
for k=0:n 

plot(k*h,x(l,k+l) , ’o') 
end; 
end; 
end ; 

plotadm=plotadm+l ; 
disp(’ ’) 

disp( ' FOLLOWING OPTIONS ARE AVAILABLE : ’) 
disp(' ’) 

disp(’ FINISH THE PROGRAM : 0 ') 

disp ( ’ DISPLAY THE TRAJECTORY, CONTROL u AND 

ACCELERATION : 1 ’) 

disp( ’ DISPLAY THE TRAJECTORY, VELOCITY, CONTROL u AND 
ACCELERATION : 2 ') 

disp( ' DISPLAY THE TRAJECTORY AND ACCELERATION ; 3 ') 
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dispC ' DISPLAY THE TRAJECTORY AND VELOCITY : 4 ') 
disp( ' DISPLAY THE TRAJECTORY AND CONTROL u : 5 ') 
dispC ' DISPLAY THE FIRST DERIVATIVE OF THE CONTROL 
u : 6 O 

disp(' DISPLAY THE CONTROLSIGNAL v : 7 ') 

disp( ' DISPLAY THE FIRST DERIVATIVE OF THE CONTROL 

w : 8 ’ ) 

disp( ' DISPLAY THE SECOND DERIVATIVE OF THE CONTROL 
w : 9 ’ ) 

spc.plot = input (' MAKE YOUR CHOICE : '); 
end; 


clear global; 
for k=2:plotadm 
delete (k) ; 
end; 
end; 


function [Q.cnt] = quad814mod(funfcn,a,b,tol) 

^Alteration of the original matlab toolbox program. 

’/.QUAD8 Numerical evaluation of an integral, higher order 
*/. method, q = QUAD8 ( ' F ' , A , B , TOL) approximates the 
*/. integral of F(X) from to B to within a relative error 

'/, of TOL. ’ F' is a string containing the name of the 

V, function. The function must return a 4*4-matrix 
*/. output value if given an input value. 
y, q = Inf is returned if an excessive recursion level 
*/, is reached indicating a possibly singular integral. 

y. qUAD8 uses an adaptive recursive Newton Cotes 8 panel 

y, rule . 

y, Cleve Moler, 5-08-88. 

*/, Copyright (c) 1984-94 by The MathWorks, Inc. 

*/, [q,cnt] = quad8(F,a,b,tol) also returns a function 

'/, evaluation count. 

7, Top level initialization, Newton-Cotes weights 
w = [3956 23552 -3712 41984 -18160 41984 -3712 23552 
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3956] / 14175 ; 

x ■ a + (0:8)*(b-a)/8; 

y, set up function call 
for i=x 

y = [y feval(funfcn,i)] ; 
end; 

*/, Adaptive, recursive Newton-Cotes 8 panel quadrature 
QO - zeros (4); 

[Q,cnt] ■ quad814stpmod(funfcn,a,b,tol,0,w,x,y,Q0) ; 

cnt = cnt + 9 ; 

end; 


function [Q,cnt] = quad814stpmod(FunFcn,a,b,tol,lev, 

w,xO,fO,qO) 

^Alteration of the original matlab toolbox program. 
XQUAD8STP Recursive function used by QUAD8. 

7. [q.cnt] = quad8stp(F,a,b,tol,lev,w,f ,q0) tries to 
y, approximate the integral of f(x) from a to b to 
7. within a relative error of tol . F is a string 
'/, containing the name of f . The remaining arguments 
y, are generated by quad8mod or by the recursion. 

*/, lev is the recursion level. 

y, w is the weights in the 8 panel Newton Cotes formula. 

y, xO is a vector of 9 equally spaced abscissa is the 

y, interval. 

*/, fO is a matrix of the 9 function values at x. 

7. qo is an approximate value of the integral. 
y, Cleve Moler, 5-08-88. 

*/, Copyright (c) 1984-94 by The MathWorks, Inc. 

LEVMAX = 10; 

'/, Evaluate function at midpoints of left and 
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right half intervals, 
x = zeros(l , 17) ; 
x(l : 2 : 17) = xO ; 

x(2 : 2 : 16) = (x0(l:8) + x0(2:9))/2; 

f C : ,1:4)- fO( : ,1 :4) ; 
for i«l:8 

f ( : ,8*i-3:8*i) = f eval(FunFcn,x(2*i) ) ; 
f ( : ,8*i+l :8*i+4) = (f0( : ,4*i+l:4*i+4)) ; 
end ; 

7 , Integrate over half intervals, 
h = (b-a)/16; 

Q1=0 ; Q2=0 ; 

for i=l:9 

Q1 * Q1 + h*w(i)*f ( : ,4*i-3:4*i) ; 
q2 = q2 + h*w(10-i)*f ( : ,69-i*4:72-i*4) ; 
end; 

q = qi + q 2 ; 

%% Recursively refine approximations, 
if norm(q - qo) > tol*norm(q) & lev <= LEVMAX 
c = (a+b)/2 ; 

[qi.cntl] = quad814stpmod(FunFcn,a,c,tol/2,lev+l, 

w,x(l:9),f(:,l:36),qi); 
[q2,cnt2] = quad814stpmod(FunFcn,c,b,tol/2,lev+l, 

w,x(9:17),f(:,33:68),q2); 

Q = Ql + Q2; 

cnt = cnt + cntl + cnt2; 

end 

end; 
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function x = mprl51knovel(spc_plot,t,n,cleargr, pointfcn, 

lambdal , lambda2 , ep) 

7 . 7 . THIS PROGRAM CALCULATES CONTROLLAWS FOR A ONE 
DIMENSIONAL TRAJECTORY •/.*/. 

VI. spc.plot DETERMINES WHICH GRAPH TO BE DISPLAYED, 
CHOSE AN INTEGER =<11 VI. 

7 . 7 . t IS THE TIME PERIOD FOR WHICH THE SYSTEM IS 
TO BE CONTROLLED VI. 

VI. n ARE THE NUMBER OF POINTS AT THE SPECIFIED 
TRAJECTORY, CHOSE t/n>l/10 7 . 7 . 

7 . 7 . cleargr CLEARS THE CURRENT WINDOW, CHOSE 1 OR 0 •/.’/. 

7 . 7 . pointfcn SPECIFIES THE TRAJECTORY *///. 

7 . 7 . lambdal AFFECTS THE INSTABILITY OF THE SYSTEM 7 . 7 . 

7 . 7 . USE lambdalOO TO AVOID NUMERICAL PROBLEMS WHEN 
DETERMINE THE MATRIX WW 7 . 7 . 

7 . 7 . Iambda2 ALSO EFFECTS THE STABILITY OF THE SYSTEM, 
CHOSE 11“=12 7 . 7 . 

7 . 7 . ep<>0 PUTS A ZERO IN THE TRANSFERFUNCTION 7 . 7 . 


global A B h t n 

7 . 7 . THE SYSTEM 7 . 7 . 

A=[ 0 1 0 0 0; 

0 lambdal 100; 

00010 ; 

00001 ; 

0000 lambda2] ; 

B= [ 0 ep 0 0 1] ' ; 

C= [ 1 0 0 0 0] ; 

7 . 7 . FUNCTION pointfcn DETERMINES THE SPECIFIED TRAJECTORY 7 . 7 . 

7 . 7 . NECESSARY FOR THE COMPOUND FUNCTION pointsinl2 7 . 7 . 
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if pointfcn= =, pointsinl2' ; 
t=5 . 2 ; 
n=26 ; 
end; 

R=f eval(pointf cn) ; 

XX CALCULATION OF THE INTEGRAL FROM 0 TO h %% 

m=48; XX NUMBER OF POINTS BETWEEN INTERPOLATION, 

CHOSE A MULTIPLE OF 6 */.’/. 

mp=m/6 ; NEEDED FOR fen spline.error THAT DETERMINES 

THE PRECISION IN THE SPLINE APPROXIMATION V/. 
tol=le-08; */,*/. THE NUMERIC ERROR TOLERANCE XX 

Mtau(: ,l:5)=zeros(5) ; 
tau=0 ; 
for j=l:m 
oldtau=tau; 
tau=oldtau+h/m; 

Mtau ( : ,5*j+l :5*j+5) = quad815inod ( ’ integrand * , oldtau , tau , tol) 

+ Mtau( : , 5*j-4:5*j) ; 
end; 

M=Mtau( : ,5*m+l :5*m+5) ; 


•/.*/. FORMING OF THE MATRICES FOR THE BLOCKDIAGONAL SYSTEM 1% 

e_Ah=expm(-A*h) ; 

Minv=inv(M) ; 

ZZ=Minv*e_Ah; 

WW=e_Ah ' *ZZ+Minv ; 


*/.*/. CONTINUOUS CONTROLLAW */.*/. 

WLD0= [ep*WW (2 , 2) +WW (5 , 2) ep*WW (2 , 3) +WW(5 ,3) 

ep*WW(2 ,4) +WW(5,4) ep*WW(2,5)+WW(5,5)] ; 
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ZLUDO- [ep*ZZ (2 , 2 ) +ZZ (5 , 2) ep*ZZ (2 , 3) +ZZ (5 , 3) 

ep*ZZ (2 ,4) +ZZ (5 , 4) ep*ZZ(2,5)+ZZ(5,5)] ; 

ZLLD0= [ep*ZZ (2 , 2) +ZZ (2 , 5) ep*ZZ (3 , 2) +ZZ (3 , 5) 

ep*ZZ(4,2)+ZZ(4,5) ep*ZZ(5,2)+ZZ(5,5)] ; 

WRD0=[ep*WW(2, 1)+WW(5, 1)] ; 

ZRUD0= [ep*ZZ (2 , 1 ) +ZZ (5 , 1 ) ] ; 

ZRLD0=[ep*ZZ(l,2)+ZZ(l,5)] ; 

VI. CONTINUOUS FIRST DIFFERENTIAL OF CONTROLLAW */.*/. 

WLDl=[ep*WW(l,2)+ep*lambdal*WW(2,2)+WW(4,2)+lambda2*WW(5,2) 
ep*WW (1,3) +ep*lambdal*WW(2 , 3) +WW (4 , 3) +lambda2*WW (5 , 3) 
ep*WW (1,4) +ep*lambdal*WW(2 , 4) +WW (4 , 4) +Iambda2*WW(5 ,4) 
ep*WW(l,5)+ep*lambdal*WW(2,5)+WW(4,5)+lajnbda2*WW(5,5)] ; 

ZLUD 1= [ep*ZZ (1,2) +ep*lambdal*ZZ (2 , 2 ) +ZZ (4 , 2) +lambda2*ZZ (5 , 2) 
ep*ZZ ( 1 , 3) +ep*lambdal*ZZ (2,3) +ZZ(4 , 3) +lambda2*ZZ (5 , 3) 
ep*ZZ ( 1 , 4) +ep*lambdal *ZZ (2,4) +ZZ (4 , 4) +lambda2*ZZ (5,4) 
ep*ZZ(l ,5)+ep*lambdal*ZZ(2,5)+ZZ(4,5)+lambda2*ZZ(5 ,5)] ; 

ZLLD 1= [ep*ZZ (2 , 1 ) +ep*lambdal *ZZ (2 , 2) +ZZ (2 , 4) +lambda2*ZZ (2 , 5) 
ep*ZZ (3 , 1 ) +ep*lambdal *ZZ (3 , 2) +ZZ (3 ,■ 4) +lambda2*ZZ (3 , 5) 
ep*ZZ (4,1) +ep*lambdal*ZZ (4 , 2) +ZZ (4 , 4) +lambda2*ZZ (4,5) 
ep*ZZ(5 , 1) +ep*lajnbdal*ZZ(5 , 2) +ZZ (5 , 4) +lambda2*ZZ (5 , 5) ] ; 

WRD1= [ep*WW(l , l)+ep*lambdal*WW(2 , 1)+WW(4, l)+lambda2*WW(5 , 1)] ; 

ZRUD1= [ep*ZZ(l , l)+ep*lambdal*ZZ(2 , 1)+ZZ(4, l)+lambda2*ZZ(5 , 1)] ; 

ZRLDl=[ep*ZZ(l,l)+ep*lambdal*ZZ(l,2)+ZZ(l,4)+lambda2*ZZ(l,5 )] ; 

VI. CONTINUOUS SECOND DIFFERENTIAL OF CONTROLLAW VI. 

WLD2= [ep*lambdal*WW ( 1 , 2) +ep*lambdal ~2*WW(2 , 2)+WW(3 , 2) + 

lambda2*WW(4, 2)+lambda2“2*WW(5 ,2) ep*lambdal*WW(l ,3) + 
ep*lambdal “ 2* WW ( 2 , 3) +WW (3 , 3) +lambda2*WW(4 , 3) + 
lambda2"2*WW(5 ,3) ep*lambdal*WW(l,4)+ep*lambdal*2* 

WW(2 ,4)+WW(3,4)+lambda2*WW(4,4)+lambda2"2*WW(5,4) 
ep*lambdal*WW ( 1 , 5) +ep*lambdal “2*WW(2 , 5) +WW (3,5)+ 
lambda2*WW(4,5)+lambda2~2*WW(5,5)] ; 

ZLUD2= [ep*lambdal*ZZ(l , 2)+ep*lambdal~2*ZZ(2,2)+ZZ(3,2)+ 

lambda2*ZZ(4,2)+lambda2“2*ZZ(5,2) ep*lambdal*ZZ(l,3)+ 
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ep*lambdal ‘ 2*ZZ (2 , 3) +ZZ (3 , 3) +lambda2*ZZ (4,3)+ 
lambda2‘2*ZZ(5,3) ep*lambdal*ZZ(l,4)+ep*lambdal‘2* 

ZZ (2 , 4) +ZZ (3 , 4) +lambda2*ZZ (4 , 4) +lambda2‘2*ZZ (5 , 4) 
ep*lambdal*ZZ(l , 5) +ep*lambdal ~2*ZZ(2 , 5)+ZZ(3, 5) + 
lambda2*ZZ(4,5)+lambda2"2*ZZ(5,5)] ; 

ZLLD2= [ep*lambdal*ZZ (2 , 1) +ep*lambdal ‘2*ZZ(2 , 2) +ZZ(2 , 3) + 

lambda2*ZZ(2,4)+lambda2‘2*ZZ(2,5) ep*lambdal*ZZ(3, 1)+ 
ep*lambdal‘2*ZZ(3,2)+ZZ(3,3)+lambda2*ZZ(3,4)+ 
lambda2‘2*ZZ(3,5) ep*lambdal*ZZ(4, l)+ep*lambdal‘2* 

ZZ (4 , 2) +ZZ (4 , 3) +lambda2*ZZ (4 , 4) +lambda2‘2*ZZ (4,5) 
ep*lambdal*ZZ(5 , l)+ep*lanibdal 2*ZZ(5,2)+ZZ(5,3)+ 
lambda2*ZZ(5 ,4)+lambda2‘2*ZZ(5,5)] ; 

WRD2= [ep*lambdal*WW(l, l)+ep*lambdal‘2*WW(2,l)+WW(3,l)+ 
lambda2*WW(4, l)+lambda2‘2*WW(5 , 1)] ; 

ZRUD2=[ep*lambdal*ZZ(l , l)+ep*lambdal‘2*ZZ(2,l)+ZZ(3, 1)+ 
lambda2*ZZ(4 , 1) +lambda2‘2*ZZ (5 , 1)] ; 

ZRLD2=[ep*lambdal*ZZ(l , l)+ep*lambdal‘2*ZZ(l,2)+ZZ(l,3)+ 
lambda2*ZZ(l,4)+lambda2‘2*ZZ(l,5)] ; 

•/,*/, CONTINUOUS THIRD DIFFERENTIAL OF CONTROLLAW */.*/. 

WLD3= [ep*lambdal ‘2*WW ( 1 , 2) + (ep*lambdal‘3+l) *WW (2 , 2)+lambda2* 
WW ( 3 , 2 ) +lambda2 * 2 *WW (4 , 2 ) +lambda2 ‘3*WW (5,2) 

ep*lambdal‘2*WW(l,3)+(ep*lambdai‘3+l)*WW(2,3)+lambda2* 

WW (3 , 3) +lambda2~ 2*WW (4 ,3) +lambda2‘3*WW (5 , 3) 

ep*lambdal‘2*WW(l,4)+(ep*lambdal‘3+l)*WW(2,4)+lambda2* 

WW (3 , 4) +lambda2‘2*WW (4,4) +lambda2‘3*WW (5,4) 

ep * 1 ambdal “ 2 * WW ( 1 , 5 ) + ( ep* lambdal “3+l)*WW(2,5)+l ambda2* 

WW (3 , 5 ) +lambda2~ 2*WW (4 , 5 ) +lambda2 ‘3*WW (5 , 5 ) ] ; 

ZLUD3= [ep*lambdal ‘2*ZZ (1,2)+ (ep*lambdal‘3+l) *ZZ (2 , 2) +lambda2* 
ZZ(3,2)+lambda2‘2*ZZ(4,2)+lambda2‘3*ZZ(5,2) 
ep*lambdal“2*ZZ(l , 3)+(ep*lambdal‘3+l)*ZZ(2,3)+lambda2* 
ZZ (3 , 3) +lambda2‘2*ZZ (4 , 3) +lambda2~3*ZZ(5 ,3) 
ep*lambdal ‘ 2*ZZ ( 1 , 4) + (ep*lambdal ‘3+1) *ZZ (2 , 4) +lambda2* 
ZZ (3 , 4) +lambda2‘2*ZZ (4 , 4) +lambda2‘3*ZZ(5 ,4) 

ep*lambdal‘2*ZZ(l,5)+(ep*lambdal‘3+l)*ZZ(2,5)+lambda2* 
ZZ(3 , 5)+lambda2‘2*ZZ(4, 5)+lambda2‘3*ZZ(5 ,5)] ; 
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ZLLD3= [ep*lambdal “2*ZZ (2 , 1) + (ep*lambdal“3+l) *ZZ(2 , 2) +lambda2* 
ZZ(2,3)+lambda2“2*ZZ(2,4)+lambda2“3*ZZ(2,5) 
ep*lambdal “2*ZZ (3,1)+ (ep*lambdal “3+1) *ZZ(3 , 2) +lambda2* 
ZZ(3,3)+lambda2“2*ZZ(3,4)+lambda2“3*ZZ(3,5) 
ep*lambdal“2*ZZ (4,1)+ (ep*lambdal "3+1) *ZZ (4,2) +lambda2* 
ZZ(4,3)+lambda2"2*ZZ(4,4)+lambda2“3*ZZ(4,5) 
ep*lambdal “2*ZZ (5 , 1) + (ep*lambdal“3+l) *ZZ(5 , 2) +lambda2* 
ZZ(5 ,3)+lambda2“2*ZZ(5,4)+lambda2“3*ZZ(5,5)] ; 
WRD3=[ep*lambdal“2*WW(l,l)+(ep*lambdal“3+l)*WW(2,l)+lambda2* 
WW ( 3 , 1 ) +lambda2 “ 2*WW (4 , 1 ) +lambda2 “3*WW (5,1)]; 

ZRUD3= [ep*lambdal“2*ZZ ( 1 , 1) + (ep*lambdal“3+l) *ZZ(2 , 1) +lambda2* 
ZZ (3 , 1 ) +lambda2 " 2*ZZ (4 , 1) +lambda2“3*ZZ(5 , 1 ) ] ; 

ZRLD3= [ep*lambdal“2*ZZ (1,1) + (ep*lambdal“3+l) *ZZ(1 , 2) +lambda2* 
ZZ(1 ,3)+lambda2"2*ZZ(l ,4)+lambda2“3*ZZ(l ,5 )] ; 

7.7. DIFFERENTIAL APPROXIMATIONS FOR THE BOUNDARY CONDITIONS 7.7. 

ydlO=(R(l , 2) -R(l , 1) )/h; 
ydll=(R(l ,3)-R(l ,2) )/b; 
ydl2=(R(l ,4)-R(l ,3) )/h; 
ydl3=(R(l ,5)-R(l ,4) )/h; 
ydln=(R(l ,n+l)-R(l ,n))/h; 
ydln_l=(R(l , n) -R(l , n-1) ) /h; 
ydln_2=(R(l ,n-l)-R(l ,n-2))/h; 
ydln_3=(R(l ,n-2)-R(l ,n-3) )/h; 
uO= (ydlO-ydl 1) /h ; 
ul=(ydll-ydl2)/h; 
u2= (ydl2-ydl3) /h ; 
un= (ydln_ 1-ydln) /h ; 
un_l=(ydln_2-ydln_l)/h; 
un_2= (ydln_3-ydln_2) /h ; 
udlO=(uO-ul)/h; 
udll=(ul-u2)/h; 
udln= (un_ 1-un) /h ; 
udln_l=(un_2-un_l)/h; 
ud20=(udl0-udll)/h; 
ud2n= (udln_ 1 -udln) /h ; 
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4/5 of the the first state vector Y,*/, 

X0= [ydlO ; uO; udlO; ud20] ; 

y,y, 4/5 of the the last state vector *1*1, 

Xn=[ydln; un; udln; ud2n] ; 

*/,*/, FORMING OF THE RIGHT HAND SIDE OF THE 
BLOCKDIAGONAL SYSTEM */,*/, 

j =1 ; 

for i=2:n 

OmegaCj , l)=-ZRLD0*R(l , i-l)+WRD0*R(l , i)-ZRUDO*R(l , i+1) ; 

j=j+i; 

OmegaCj , 1) =-ZRLDl*R( 1 , i-l) +WRD1*R(1 , i)-ZRUDl*R(l , i+1) ; 

j=j +1 ; 

OmegaCj , 1 ) =-ZRLD2*R (1 , i-l) +WRD2*R C 1 , i ) -ZRUD2*RC 1 , i+ 1 ) ; 

j=j +1 ; 

OmegaCj ,l)=-ZRLD3*RCl,i-l)+WRD3*RCl,i)-ZRUD3*RCl, i+1) ; 
end; 

Omega C 1 , 1 ) =Omega C 1 , 1 ) -ZLLDO *X0 ; 

0megaC2, l)=0megaC2, 1)-ZLLD1*X0; 

Omega C3, l)=0megaC3, 1)-ZLLD2*X0; 

0megaC4, l)=0megaC4, 1)-ZLLD3*X0; 

0megaC4*n-7, l)=0mega(4*n - 7 > l)-ZLUD0*Xn; 

Omega C4*n-6 , l)=0megaC4*n-6, l)-ZLUDl*Xn; 

Omega C4*n-5 , l)=0megaC4*n-5 , l)-ZLUD2*Xn; 

Omega C4*n-4 , 1 ) =0mega C4*n-4 , 1 ) -ZLUD3*Xn ; 

•/.*/. FORMING OF THE LEFT HAND SIDE OF THE 
BLOCKDIAGONAL SYSTEM */,*/, 

for i=l:4 

DDCl,i)=-WLDOCl,i) ; 

DDC2,i)— WLDl(l.i) > 

DD(3,i)=-WLD2Cl,i) ; 

DDC4,i)=-WLD3Cl,i) ; 

if n>2 
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DD ( 1 , i +4) =ZLUDO ( 1 , i ) ; 

DD(2 , i+4)=ZLUDl (1 , i) ; 

DD (3 , i+4) =ZLUD2 ( 1 , i ) ; 

DD (4 , i+4) =ZLUD3 ( 1 , i) ; 
end; 
end; 

for i=2:n-2 
for j=~l : 2 

DD (4*i-3,4*i-6+j ) =ZLLD0 ( 1 , j +2) ; 

DD (4*i-3 , 4*i-2+j ) =-WLD0 ( 1 , j +2) ; 

DD (4*i-3 , 4*i+2+ j ) =ZLUD0 ( 1 , j +2) ; 
DD(4*i-2,4*i-6+j)=ZLLDl(l, j+2) ; 
DD(4*i-2,4*i-2+j)=-WLDl(l, j+2) ; 

DD (4*i-2 ,4*i+2+ j ) =ZLUD1 (l,j+2); 
DD(4*i-l,4*i-6+j)=ZLLD2(l, j+2) ; 

DD (4*i-l ,4*i-2+j ) =-WLD2 ( 1 , j +2) ; 
DD(4*i-l,4*i+2+j)=ZLUD2(l, j+2) ; 

DD(4*i , 4*i-6+j )=ZLLD3 (1 , j +2) ; 

DD(4*i , 4*i-2+j)=-WLD3(l, j+2) ; 

DD(4*i , 4*i+2+j )=ZLUD3 (1 , j +2) ; 

end; 
end; 
if n>2 

for i=l :4 

DD(4*n-7,4*n-12+i)=ZLLD0(l ,i) ; 
DD(4*n-6,4*n-12+i)=ZLLDl(l ,i) ; 
DD(4*n-5,4*n-12+i)=ZLLD2(l,i) ; 
DD(4*n-4,4*n-12+i)=ZLLD3(l,i) ; 

DD (4*n-7 ,4*n-8+i)=-WLD0 ( 1 , i) ; 

DD (4*n-6 ,4*n-8+i)=-WLDl ( 1 , i) ; 

DD(4*n-5 ,4*n-8+i)=-WLD2(l , i) ; 

DD (4*n-4,4*n-8+ i ) =-WLD3 ( 1 , i ) ; 
end; 
end; 

DD=sparse(DD) ; ’/.•/. SQUEEZING OUT ALL ZERO ELEMENTS 

FROM MATRIX DD XX 
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%•/. GAUSSELIMINATION TO PRODUCE AN UPPER TRIANGULAR SYSTEM 

for k=2:4:4*n-10 
for l=k:k+6 

if DD(k-l ,k-l)"=0 

zd=DD(l,k-l)/DD(k-l,k-l) ; 

DD(1, : )=DD(1, : )-zd*DD(k-l , : ) ; 

Omega(l, l)=0mega(l, l)-zd*Omega(k-l,l) ; 
end; 
end; 

for l=k+l:k+6 
if DD(k,k) “=0 

zd=DD(l,k)/DD(k,k) ; 

DD(1, :)*DD(1, :)-zd*DD(k, :); 

Omega(l , l)=0mega(l , 1) -zd*0mega(k, 1) ; 
end; 
end; 

for l=k+2:k+6 

if DD(k+l,k+l)~=0 

zd=DD(l,k+l)/DD(k+l,k+l) ; 

DD(1, : ) =DD (1 , : ) -zd*DD (k+ 1 , : ) ; 

Omega (1 , 1) =0mega(l , 1) -zd*Omega(k+l , 1) ; 
end; 
end; 

for l=k+3:k+6 

if DD(k+2,k+2)"=0 

zd=DD(l,k+2)/DD(k+2,k+2) ; 

DD (1 , : ) =DD (1 , : ) -zd*DD (k+2 , : ) ; 

Omega (1 , 1) =0mega(l , 1) -zd*0mega(k+2 , 1) ; 
end; 
end; 
end; 

k=4*n-7 ; 
if DD(k,k ) ~-0 

zd=DD(k+l ,k)/DD(k,k) ; 

DD(k+l , : ) =DD(k+l , : )-zd*DD(k, : ) ; 
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Omega (k+ 1,1) =Omega (k+ 1 , 1 ) -zd*Omega (k , 1) ; 
zd=DD(k+2,k)/DD(k,k) ; 

DD (k+2 , : ) =DD (k+2 , : ) -zd*DD (k , : ) ; 

0mega(k+2 , 1) =0mega(k+2 , 1) -zd*Omega(k, 1) ; 
zd=DD(k+3,k)/DD(k,k) ; 

DD (k+3 , : ) =DD (k+3 , : ) -zd*DD (k , : ) ; 

Omega (k+3 , l)=0mega(k+3, l)-zd*Omega(k, 1) ; 
end; 

if DD(k+l,k+l)"=0 

zd=DD(k+2,k+l)/DD(k+l,k+l) ; 

DD (k+2 , : ) =DD (k+2 , : ) -zd*DD (k+ 1 , : ) ; 

0mega(k+2 , 1) =0mega(k+2 , 1) -zd*Omega(k+l , 1) ; 
zd=DD (k+3 , k+ 1 ) /DD (k+ 1 ,k+l) ; 

DD (k+3 , : )=DD (k+3 , : ) -zd*DD (k+1 , : ) ; 

Omega (k+3 , l)=0mega(k+3, 1) -zd*Omega(k+l , 1) ; 
end; 

if DD(k+2 ,k+2) "=0 

zd=DD (k+3 , k+2) /DD (k+2 , k+2) ; 

DD (k+3 , : ) =DD (k+3 , : ) -zd*DD (k+2 , : ) ; 

Omega (k+3 , 1) =0mega(k+3 , 1) -zd*0mega(k+2 , 1) ; 
end; 

*/,•/, BACKSUBSTITUTION TO SOLVE FOR THE STATEVECTORS */.*/. 

ud2 (n- 1 ) =Omega(4*n-4 , 1 ) /DD (4*n-4 , 4*n-4) ; 

udl (n- 1 ) = (0mega(4*n-5 , 1) -DD (4*n-5 ,4*n-4) *ud2 (n-1) ) / 

DD(4*n-5 ,4*n-5) ; 

u (n - 1 ) = (0mega(4*n~6 ,1) - DD (4*n~6 ,4*n~4) *ud2 (n~l) ~ 

DD (4*n-6 , 4*n-5) *udl (n- 1 ) ) /DD (4*n-6 , 4*n-6) ; 

ydl (n- 1 ) = (Omega(4*n-7 , 1 ) -DD (4*n-7 , 4*n-4) *ud2 (n-1)- 

DD (4*n-7 , 4*n-5 ) *udl (n- 1 ) -DD (4*n-7 , 4*n-6) *u(n-l) ) / 

DD(4*n-7 ,4*n-7) ; 

if n>2 

for k=n-2 : -1 : 1 

ud2 (k) = (Omega (4*k , 1 ) -DD (4*k ,4*k+4) *ud2 (k+1) - 
DD(4*k,4*k+3)*udl(k+l)-DD(4*k,4*k+2)*u(k+l)- 
DD (4*k , 4*k+l) *yd 1 (k+1 ) ) /DD (4*k , 4*k) ; 
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udl (k) = (0mega(4*k- 1 , 1 ) ■ -DD (4*k-l , 4*k+4) *ud2 (k+ 1) - 
DD (4*k- 1 , 4*k+3) *ud 1 (k+ 1 ) -DD (4*k- 1 , 4*k+2) *u (k+1 ) - 
DD (4*k- 1 , 4*k+ 1) *ydl (k+1 ) -DD (4*k- 1 , 4*k) *ud2 (k) ) / 
DD(4*k-l ,4*k-l) ; 

u (k) = (0mega(4*k-2 , 1 ) -DD (4*k-2 ,4*k+4) *ud2 (k+1) - 
DD (4*k-2 ,4*k+3) *udl (k+1) -DD(4*k-2 ,4*k+2) *u(k+l) - 
DD (4*k-2 , 4*k+ 1 ) *ydl (k+1 ) -DD (4*k-2 , 4*k) *ud2 (k) - 
DD(4*k-2 ,4*k-l)*udl (k) ) /DD(4*k-2 ,4*k-2) ; 
yd 1 (k) = (Omega (4*k-3 , 1 ) -DD (4*k-3 , 4*k+4) *ud2 (k+ 1 ) - 
DD (4*k-3 , 4*k+3) *udl (k+ 1 ) -DD (4*k-3 , 4*k+2) *u (k+ 1 ) - 
DD (4*k-3 , 4*k+ 1 ) *ydl (k+ 1 ) -DD (4*k-3 , 4*k) *ud2 (k) - 
DD(4*k-3 ,4*k-l)*udl (k)-DD(4*k-3 ,4*k-2)*u(k))/ 

DD (4*k-3 , 4*k-3) ; 
end; 
end; 


*/.y, MAKING OF THE STATEVECTORS 

x(: ,1)-CR(1.1); X03 ; 
x( : ,n+l)=[R(l,n + l) ; Xn] ; 
for i=l:n-l 

x(: ,i+l)=[R(l»i + l) ; ydl(i); u(i); udl(i); ud2(i)J ; 
end; 

if cleargr 
elf ; 

hold on; 
grid on; 
end; 

•/,*/. PLOTTING OF THE CALC/SPEC TRAJECTORY, VELOCITY, 
CONTROLSIGNAL AND ACCELERATION 

plotadm=0 ; 
while spc_plot 
if plotadm 

if spc_plot<12 
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figure; 
hold on; 
grid on; 

title(['0ne dim. traj : mpr!51knovel .m 11 = 

' ,num2str(lambdal) , ' 12 = ' ,num2str(lambda2) , ' ep = 

' f num2str(ep) , ' n = ' ,num2str(n)] ) 

xlabel( ['Total error : ’ ,num2str(Total_error) , 

' Point error : ' ,num2str(Point_error) , ' Graph : 
' ,num2str(spc_plot) , ' Time t']) 
if spc_plot<6 
for k=0:n 

plot(k*h,x(l ,k+l) , 'o') 
end; 
end; 

end; 

end; 

break adm=0 ; 

fadm=0; */,•/. NORMALLY fadm=0, NECESSARY FOR WRITING OF 
TEXTS IN THE PLOT •/.*/, 

for j=0:m 

eAtau=expm(A*j*h/m) ; 

for i=f adm:n-l •/.*/. CHOOSE e.g 2:n-3 TO AVOID 

PROBLEMS AT THE ENDPOINTS 7.’/. 
entry ( : , i+l)=eAtau*(x( : , i+l)+Mtau( : ,5*j+l :5*j+5)* 
Minv*(e_Ah*x( : , i+2) -x( : , i+1)) ) ; 
csignvec( : , i+l)=B'*expm(-A' *j *h/m)*Minv* 

(e_Ah*x( : , i+2)-x( : , i+1)) ; 
if j==0 

if i==f adm ; 

entryl=entry(l,fadm+l) ; 
entry2=entry (2,f adm+1) ; 
entry3=entry(3,fadm+l) ; 
entry4=entry(4,fadm+l) ; 
entry5=entry(5,fadm+l) ; 
csignvecl=csignvec(l , 1) ; 
end;0 
end; 
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if rem(j,mp)==0 
if j <=m-mp 

t ra j ect ( 1 , j /mp*n+ ( i+ 1 ) ) =entry ( 1 , i+ 1) ; 
end; 
end; 
if j==m 
if i==n-l 

traj ect (1 , 6*n+l)=entry(l, i+1) ; 
end; 
end; 

if spc_plot==l 

plot (i*h+j*h/m, entry (1 , i+1) ) VI. TRAJECTORY 

plot(i*h+j*h/m,entry(3,i+l) , ’ . ’) */•'/• CONTROL u '/,*/. 

plot (i*h+j*h/m, lambdal*entry(2, i+l)+entry (3, i+l)+ 
ep*cs ignvec (1 , i+1 ) , ' . » ) ’/.*/. ACCELERATION •/.'/. 

elseif spc_plot==2 

plot (i*h+j*h/m, entry(l , i+1) ,’•' ) VI. TRAJECTORY 

plot (i*h+j*h/m, entry (2, i+1) , ’ . ') VI. VELOCITY VI. 

plot(i*h+j*h/m,entry(3,i+l) , ' . ') VI. CONTROL u VI. 
plot (i*h+j*h/m, lambdal*entry(2, i+l)+entry(3, i+l)+ 
ep*csignvec(l , i+1) , ' . ' ) VI. ACCELERATION '1.1, 
elseif spc_plot==3 

plot(i*h+j*h/m,entry(l,i+l) , ’ . ’) VI. TRAJECTORY VI, 
plot (i*h+j*h/m,lambdal*entry(2 , i+l)+entry(3, i+l)+ 
ep*csignvec(l , i+1) , ' . ' ) VI, ACCELERATION 1,1, 
elseif spc_plot==4 

plot(i*h+j*h/m,entry(l,i+l) , ' • ') VI, TRAJECTORY VI, 
plot (i*h+j*h/m , entry (2 , i+1) ) VI. VELOCITY 1,1, 

elseif spc_plot==5 

plot (i*h+j*h/m, entry (1 , i+1) ) VI, TRAJECTORY VI, 

plot(i*h+j*h/m,entry(3,i+l) , ' . ’) VI. CONTROL u '/,/, 
elseif spc_plot==6 

VI. CONTROL DER u-dot VI. 
plot (i*h+j*h/m, entry (4, i+1) ,’ . ’) 
elseif spc_plot==7 

VI. CONTROL DERx2 u-dot-dot VI. 
plot (i*h+j*h/m , entry (5 , i+1) , ' . ’) 
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elseif spc_plot==8 

csignvec( : , i+l)=B'*expm(-A'*j*h/in)*Minv* 
(e_Ah*x( : ,i+2)-x( : , i+1) ) ; 

*/.*/. CONTROLS I GNAL w 7 . 7 . 
plot (i*h+j *h/m, csignvec (1 , i+1) , ' . ') 
elseif spc_plot==9 

csdlvec( : ,i+l)=B , *A'*expin(-A , *j*h/m)*Minv* 
(e_Ah*x( : ,i+2)-x( : ,i+l)) ; 

•/.*/. CONTROL DER w-dot VI. 
plot(i*h+j*h/m,csdlvec(l,i+l) , ’ . ') 
if j==0 

csdvecl=csdlvec(l , 1) ; 
end; 

elseif spc_plot==10 

csd2vec( : ,i+l)=B'*A , *A’*expm(-A , *j*h/m)*Minv* 
(e_Ah*x( : ,i+2)-x( : ,i+l)) ; 

•/.*/. CONTROL DERx2 w-dot-dot 7 , 7 . 
plot(i*h+j*h/m,csd2vec(l,i+l) , ’ . ') 
if j==0 

csdvec2=csd2vec(l , 1) ; 
end; 

elseif spc_plot==l 1 

csd3vec( : , i+l)=B ’ *A’ *A ' *A ' *expm(-A , *j*h/m)* 
Minv*(e_Ah*x( : ,i+2)-x(: ,i+l)) ; 

7 , 7 . CONTROL DERx3 w-dot-dot-dot 7 . 7 . 
plot(i*h+j*h/m,csd3vec(l,i+l) , ' . 
if j”0 

csdvec3=csd3vec(l , 1) ; 
end; 
else 

disp(' ’) 

disp( ' NOT A VALID CHOICE ’) 
breakadm=l ; 
end; 

if breakadm 
break; 
end; 
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end; 

if breakadm 
break; 
end; 
end; 

if spc_plot<12 
ax = ax is ; 
ax2=ax(l ,2) ; 
if ax2<l . 5 

xpos=3/4*0.2; 
elseif ax2>=5 
xpos=3/4; 
else 

xpos=0 . 25 ; 
end; 

if spc_plot<3 
ax4=ax(l ,4) ; 
if ax4<l 

ax4=ax4*10; 
if ax4<l 

ax4=ax4*10; 

end; 

end; 

if ax4>10 
ax4=ax4/10; 
if ax4>10 
ax4=ax4/ 10 ; 
end ; 
end; 

if rem(ax4,2)==0 

yadd= (ax (1 ,4)/4)*l/5; 
else 

yadd= (ajc(l ,4)/3)*l/5; 
end; 

ypos=lambdal*entry2+entry3+ep*csignvecl 

text(-xpos ,ypos, ['ACC']) ; 
if abs(ypos-entry3)>yadd 
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text (-xpos, entry3, ['CS u']); 
elseif ypos-entry3>0 

text (-xpos, entry3-yadd, ['CS u'] ) ; 
else 

text(-xpos,entry3+yadd, [’CS u']) ; 
end; 

if spc_plot==2 

text(-xpos,entry2, C'VEL’]) ; 
end; 
end; 

if spc_plot==3 

ypos=lambdal*entry2+entry3+ep*csignvecl; 

text (-xpos, ypos, ['ACC']) ; 
end; 

if spc_plot==4 

text (-xpos,entry2, [ J VEL'] ) ; 
end; 

if spc_plot==5 

text(-xpos,entry3, ['CS u’]); 
end; 

if spc_plot==6 

text (-xpos, entry4, ['udl'] ) ; 
end; 

if spc_plot==7 

text (-xpos, entry5, ['ud2']) ; 
end; 

if spc_plot==8 

text (-xpos, csignvecl, [’CS w']); 
end ; 

if spc_plot==9 

text (-xpos , csdvecl , [' wdl '] ) ; 
end; 

if spc_plot==10 

text (-xpos , csdvec2 , [ ' wd2 ’ ] ) ; 
end; 

if spc_plot==ll 

t ext ( -xpos, csdvec3, ['wd3']) ; 
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end; 

end; 


•/.*/. ESTIMATION OF THE ACCURACY IN THE 
SPLINE APPROXIMATION */.*/. 

if plotadm==0 

[Total_error , Average_error , Point .error ,End_point .error] - 

spline_error(pointf cn,R,traject) ; 

Total.error 

Average.error 

Point.error 

End.point.error 

title( ['One dim. traj : mprl51knovel .m 11 = 

' , num2str (lambdal) , ’ 12 = ' ,num2str(lambda2) ,'ep- 

» ,num2str(ep) , ’ n = * ,num2str(n)] ) 

xlabel( [’Total error : ' ,num2str (Total.error) , 

> Point error : ’ ,num2str(Point_error) , ' Graph : 

’ ,num2str(spc_plot) , ' Time t']) 
if spc_plot<6 
for k=0:n 

plot(k*h,x(l,k+l) , ’o’) 
end; 
end; 
end; 

plot adm=pl ot adm+ 1 ; 
disp(' ’) 

disp( ’ FOLLOWING OPTIONS ARE AVAILABLE : ') 
disp ( ' ') 

disp( ' FINISH THE PROGRAM : 0 ’) 

disp( ' DISPLAY THE TRAJECTORY, CONTROL u AND 

ACCELERATION : 1 ') 

disp( ’ DISPLAY THE TRAJECTORY, VELOCITY, CONTROL u AND 
ACCELERATION : 2 ’) 

disp( ’ DISPLAY THE TRAJECTORY AND ACCELERATION : 3 ') 
disp( ' DISPLAY THE TRAJECTORY AND VELOCITY : 4 ') 
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disp( ' DISPLAY THE TRAJECTORY AND CONTROL u : 5 ’) 

disp( ' DISPLAY THE FIRST DERIVATIVE OF THE CONTROL u : 6 ') 

disp( ' DISPLAY THE SECOND DERIVATIVE OF THE 

CONTROL u : 7 ’ ) 

disp( ’ DISPLAY THE CONTROLSIGNAL w : 8 ') 

disp( ' DISPLAY THE FIRST DERIVATIVE OF THE CONTROL w : 9 ') 
disp( ' DISPLAY THE SECOND DERIVATIVE OF THE 
CONTROL w : 10 ') 

disp( ’ DISPLAY THE THIRD DERIVATIVE OF THE 
CONTROL w : 11 ') 

spc.plot = input ( ’ MAKE YOUR CHOICE : ’); 
end ; 


clear global; 
for k=2:plotadm 
delete (k) ; 
end; 
end; 


function [Q,cnt] = quad815mod(funfcn,a,b,tol) 

’/.Alteration of the original matlab toolbox program. 

7.QUAD8 Numerical evaluation of an integral, higher order 
*/. method, q = QUAD8 ( ' F ' , A , B , TOL) approximates the 
*/. integral of F(X) from to B to within a relative error 

*/. of TOL. *F' is a string containing the name of the 

7 function. The function must return a 5*5-matrix 
*/, output value if given am input value. 

*/, Q = Inf is returned if an excessive recursion level 
*/. is reached indicating a possibly singular integral. 

*/, qUAD8 uses an adaptive recursive Newton Cotes 8 panel 

y. rule . 

y. Cleve Moler , 5-08-88. 

•/. Copyright (c) 1984-94 by The MathWorks, Inc. 

•/. [q,cnt] = quad8(F,a,b,tol) also returns a function 

*/, evaluation count. 

*/, Top level initialization, Newton-Cotes weights 
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w = [3956 23552 -3712 41984 -18160 41984 -3712 23552 
39561/14175; 

x - a + (0 :8)*(b-a)/8; 

X set up function call 
for i=x 

y - [y f eval (funf cn, i)] ; 
end; 

X Adaptive, recursive Newton-Cotes 8 panel quadrature 
Q0 = zeros(5) ; 

[Q,cnt] = quad815stpmod(funfcn,a,b,tol,0,w,x,y,q0) ; 

cnt = cnt + 9 ; 

end; 


function [Q,cnt] = quad815stpmod(FunFcn,a,b,tol,lev, 

w,xO,fO,QO) 

XAlteration of the original matlab toolbox program. 
XQUAD8STP Recursive function used by QUAD8. 

*/. [q.cnt] » quad8stp(F,a,b,tol,lev,w,f ,q0) tries to 
X approximate the integral of f(x) from a to b to 
X within a relative error of tol. F is a string 
X containing the name of f . The remaining arguments 
X are generated by quad8mod or by the recursion. 

X lev is the recursion level. 

X w is the weights in the 8 panel Newton Cotes formula. 

X xO is a vector of 9 equally spaced abscissa is the 

X interval. 

X fO is a matrix of the 9 function values at x. 

X qo is an approximate value of the integral. 

X Cleve Moler, 5-08-88. 

X Copyright (c) 1984-94 by The MathWorks, Inc. 

LEVMAX = 10; 

X Evaluate function at midpoints of left and 



10 MATLAB PROGRAMS 


102 


right half intervals, 
x ■ zeros(l , 17) ; 
x(l :2 : 17) = xO; 

x(2:2: 16) = (x0(l:8) + x0(2:9))/2; 


f ( : , 1 :5)“ f0(: ,1:5) ; 
for i=l:8 

f ( : ,10*i-4: 10*i) = feval(FunFcn,x(2*i)) ; 
f (: ,10*i+l:10*i+5) = (f 0( : ,5*i+l : 5*i+5)) ; 
end; 

V, Integrate over half intervals, 
h = (b-a)/16; 

Q1=0;Q2=0; 
for i=l:9 

qi = qi + h*w(i) *f ( : , 5*i-4 : 5*i) ; 
q2 = q2 + h*w(10-i)*f ( : ,86-i*5 : 90-i*5) ; 
end; 

q = qi + q 2 ; 

*/,*/, Recursively refine approximations, 
if norm(q - qo) > tol*norm(q) & lev <= LEVMAX 
c = (a+b)/2; 

[qi.cntl] = quad815stpmod(FunFcn,a,c,tol/2,lev+l, 

w,x(l :9) ,f (: ,1:45) ,qi) ; 

[q2 ,cnt2] = quad815stpmod(FunFcn,c,b,tol/2,lev+l, 

w,x(9:17),f(:, 41:85), q2); 

q = Qi + q 2 ; 

cnt = cnt + cntl + cnt2; 

end 

end; 
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function res = integrand(v) 

*/,*/, THIS SUBPROGRAM DETERMINE THE INTEGRAND OF THE MATRIX 
CONSTANT M XX 

global A B 

e_AvB=expm(-A*v)*B; 
res = e_AvB*e_AvB' ; 
end; 


function [Total_error,Average_error .Point .error , 

End.point. error] =spline_error(pointfcn,R, traj ect) 

global h t n 

Total.error=0; 

Point_error=0; 

n=6*n; 

Rp=f eval(pointfcn) ; 
n=n/6 ; 
h=t/n; 
for i=0 :n-l 
for j=0:5 

Total_error=Tot al.error+abs (traj ect (1 , j *n+l+i) - 

Rp(l , i*6+j+l) ) ; 
if j==0 

Po int_error=Point_error+abs (traj ect (l,i+l)-R(l , i+l) ) ; 

end; 

end; 

end; 

Average_error=Total_error/ (6*n) ; 

End.point _error=abs (traj ect (1 , 6*n+l) -R(l ,n+l) ) ; 
end; 
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function R=pointexpll 
*/, R * l*(n+l) -matrix, 
global h t n 

•a TIME interval for renewal of the trajectory 7.*/. 

h=t/n; 

for j=0:h:n*h 

R(1 , j/h.+l) = ( 1-exp (-3/2* (j-t/2) ) )/ (l+exp(-3/2*(j-t/2)) ) ; 
end; 
end; 


function R=pointsinl2 
% R = l*(n+l) -matrix, 
global h t n 

y.7. TIME INTERVAL FOR RENEWAL OF THE TRAJECTORY V/. 
h=t/n; 

y.7, I APOLOGIZE FOR THE "SMART" PROGRAMING %% 
adm=n*3/26; 
for j = 1 : a dm 
R(l,j)=0; 
end; 

for j=0 :h: (n-2*adm)*h 

R(1 , j/h+adm+l)=l/10*(l+sin(-pi/2+j*pi/2)) ; 
end; 

for j = 1 :adm 

R(1 , j-adm+n+l)=0; 
end; 
end; 
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function R=pointstepl 
*!, R - l*(n+l) -matrix, 
global h t n 
h=t/n; 

for j=0:h:n*h 

R(l, j/h+l)=l/2*(l+sign(j-(t/2+0 .01))) ; 
end; 
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